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SECTION  1 


INTRODUCTION  AND  SUMMARY 


1 . 1  BACKGROUND . 


During  the  Second  World  War,  both  sides  quickly  discovered 
that  fire,  rather  than  blast,  was  the  most  efficient  means  of 
damaging  large-area  targets  such  as  cities.  In  the  latter  stages 
of  the  strategic  air  offensive  against  Germany,  for  example, 
approximately  70%  of  the  total  tonnage  of  bombs  dropped  were 
incendiaries  [1].  In  certain  instances  where  a  large  number  of 
incendiaries  were  successfully  dropped  within  a  relatively  small 
area,  the  resulting  fires  merged  into  one,  termed  a  firestorm. 
Firestorms  are  characterized  by  the  nearly  complete  burning  of 
combustibles  over  a  limited  area,  rather  than  partial  combustion 
over  a  wide  front.  They  generate  severe  surface  effects  due  to 
high  temperatures  and  winds,  and  are  associated  with  the  lofting 
of  soot  and  debris  to  extreme  altitudes. 

A  nuclear  airburst  occurring  below  an  altitude  of  100,000 
feet  delivers  approximately  35%  of  its  energy  yield  in  the  form 
of  electromagnetic  radiation  at  infrared  and  visible 
wavelengths  [2].  Since  the  incendiary  effect  of  this  radiation 
is  concentrated  within  a  compact  area,  nuclear  bursts  could  serve 
as  efficient  mechanisms  for  initiating  firestorms.  While  the 
burning  radii  for  the  fires  following  the  Hiroshima  and  Nagasaki 
attacks  were  of  order  1  km  and  those  for  the  largest  firestorms 
occurring  in  the  European  theater  ranged  up  to  3  km,  the  expected 
burn  radius  for  a  nominal  1-Mt  airburst  is  10  km.  Due  to  the 
unprecedented  size  of  the  postulated  fires  and  the  limitations  of 
simplified  analytical  plume  models,  numerical  simulations  are 
necessary  to  provide  reliable  estimates  of  the  near-surface 
environments  and  particulate  clouds  generated  by  mass  fires 
following  a  nuclear  airburst. 


1.2  APPROACH. 


In  this  study,  several  simulations  were  made  with  the  DICE 
code  in  order  to  investigate  the  dependence  of  the  predicted 
firestorm  environments  on  a  variety  of  physical  and  computational 
factors.  DICE  is  an  axisymmetric  two-dimensional  numerical  model 
based  on  the  compressible  flow  equations  [3].  The  grid  spacing 
can  vary  both  horizontally  and  vertically,  with  each  cell 
representing  an  average  over  an  annular  region  (see  Fig.  1).  The 
model  predicts  the  flow  of  air  obeying  the  AFWL  equation  of 
state,  and  can  include  solid,  liquid,  or  gaseous  contaminants 
(e.g.,  soot  and  water  species),  which  exchange  heat  and  momentum 
with  the  air.  The  process  of  combustion  is  modeled  by  adding 
heat  at  a  prescribed  rate  to  the  lowest  layer  of  computational 
cells  over  the  burning  area.  The  initial  atmospheric  temperature 
profile  was  taken  from  the  U.S.  Standard  Atmosphere,  with  the 
tropopause  located  at  an  altitude  of  11  km  (see  Fig.  2). 

Section  2  contains  a  brief  overview  of  the  model.  In 
Section  3,  an  inviscid  version  of  the  model  is  used  to 
investigate  the  effects  of  variations  in  the  computational  zone 
size,  burning  rate,  burning  area,  and  multi-phase  physics.  In 
Section  4  a  simplified  turbulence  model  is  introduced,  and  the 
influence  of  turbulence  on  the  temperatures  and  stabilization 
heights  of  thermal  plumes  is  discussed.  Section  5  describes  the 
validation  of  the  model  by  the  successful  simulation  of  a  test 
fire  conducted  at  the  Meteotron  research  facility  in  France.  In 
Section  6  the  question  of  swirl  is  addressed,  and  it  is  shown 
that  significant  rotational  velocities  can  be  generated  by  large- 
area  fires. 


Multiphase  DICE  Code  Numerical  Simulation  Approach 
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Figure  I.  Schematic  illustration  of  DICE  grid  and  cell  geometry. 
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1 . 3  SUMMARY . 


Inviscid  Cases 


An  inviscid  version  of  the  model  was  used  to  study  the 
effects  of  variations  in  the  burn  rate,  burn  area,  computational 
zone  size,  and  multi-phase  physics  on  the  predicted  firestorm 
environments  (see  Table  1).  Two  different  time-dependent  burning 
rates  were  used  (see  Fig.  3)s  a  "fast"  scenario  in  which  the 
maximum  burn  rate  of  250  kw/m2  (representative  of  firestorms) 
was  maintained  for  30  minutes,  and  a  "slow"  scenario  in  which  the 
maximum  burn  rate  of  62.5  kw/m2  (representative  of  group  fires) 
was  maintained  for  120  minutes. 

The  main  conclusions  drawn  from  the  inviscid  cases  were  as 
follows : 

(1)  When  the  lower  layer  (where  the  combustion  heating  is  added) 
is  given  a  depth  of  only  100  m  and  the  "fast"  burning  scenario  is 
used  (Cases  301  and  502),  temperatures  in  excess  of  1000°K  are 
generated  (see  Fig.  4,  which  shows  a  time  history  of  the  maximum 
temperature  for  Case  502).  The  large  buoyancy  which  results 
causes  intense,  jet-like  plumes  to  rise  along  the  axis  to 
altitudes  in  excess  of  40  km,  where  they  become  over-dense  and 
fall  back  downwards.  This  tendency  of  the  plumes  to  overshoot 
their  stabilization  heights  leads  to  solutions  which  are  quasi- 
periodic  (see  Fig.  4;  see  also  Figs.  5a-b,  which  show  time 
histories  of  the  inflow  velocities  at  various  locations  in  the 
lower  layers  of  Cases  301  and  502),  and  causes  the  formation  of 
two  distinct  regions  of  lofted  soot  (see  Figs.  6a-b,  which  show 
the  positions  reached  after  one  hour  by  massless  tracers  released 
every  five  minutes  from  within  the  burning  area  of  Cases  301 
and  502),  one  in  the  stratosphere  near  the  axis,  and  another 
extending  radially  outwards  near  the  altitude  of  the  tropopause. 


Table  1 


Inviscid  firestorm  simulation  cases. 


CASE 

BURN 

RATE 

BURN 

RADIUS  (KM) 

MINIMUM 
AZ  (km) 

MINIMUM 
Ar  (km) 

SOOT 
(%  Of 

WATER  VAPOR 
mass  burned) 

301 

fast 

10 

0.1 

0.33 

0 

0 

502 

fast 

10 

0.1 

0.67 

0 

0 

504 

slow 

10 

0.1 

0.67 

0 

0 

601 

fast 

10 

1.0 

1.0 

0 

0 

603 

fast 

10 

1.0 

1.0 

5 

50 

604 

slow 

10 

1.0 

1.0 

5 

50 

605 

fast 

30 

1.0 

1.0 

5 

50 

606 

fast 

3 

1.0 

1.0 

5 

50 

608 

fast 

10 

1.0 

1.0 

5 

50  + 

ambient 


re  4.  Time  histories  of  maximum  temperatures  for  Cases  502,  504,  and  601 


Radial  velocities  in  lowest  layer  of  Figure  5b.  Radial  velocities  in  lowest  layer  of 
computational  cells  at  various  radial  computational  cells  at  various  radial 
distances  for  Case  301.  distances  for  Case  502. 


(2)  When  the  lower  layer  is  given  a  depth  of  1  km  (the  600- 
series  in  Table  1),  the  combustion  heating  is  effectively  "mixed'' 
over  the  lowest  1  km  of  the  domain,  leading  to  much  cooler 
temperatures  (see  the  plot  in  Fig.  4  of  the  maximum  temperatures 
for  Case  601,  which  differed  from  Case  502  only  in  the  changed 
zoning) .  The  reduced  buoyancy  lessens  the  tendency  of  the  plume 
to  oversnoot,  and  it  approaches  a  steady-state  with  a  maximum 
altitude  of  about  18  km  by  30  minutes  into  the  run  (see  Figs. 
7a-b,  which  show  velocity  vector  plots  of  Case  601  at 

30  and  45  minutes).  As  a  result,  the  upper  cloud  of  soot  tracers 
which  developed  in  Cases  301  and  502  does  not  appear  in  the 
tracer  plot  for  Case  601  (see  Fig.  6d) ,  although  the  lower  tracer 
cloud  extending  outwards  along  the  tropopause  is  still  present. 

(3)  Additional  information  is  gained  by  using  the  multiphase 
capabilities  of  the  model  to  predict  concentrations  of  soot, 
water  vapor,  liquid  water,  and  ice  within  the  grid.  For 
Cases  603-608  in  Table  1,  5%  of  the  fuel  burned  is  converted  to 
soot,  and  50%  to  water  vapor;  Cases  603-606  neglect  ambient  water 
vapor,  while  Case  608  assumes  an  initial  relative  humidity  which 
varies  linearly  from  60%  at  the  surface  to  zero  at  an  altitude  of 
30  km.  Case  603,  using  the  same  parameters  as  the  single-phase 
Case  601,  predicts  a  soot  concentration  field  at  one  hour  (see 
Fig.  8)  which  resembles  the  tracer  cloud  from  Case  601  (see 

Fig.  6d),  but  also  shows  that  the  highest  soot  concentrations  are 
located  near  the  axis,  and  that  considerable  vertical  dispersion 
occurs  within  the  cloud.  The  ice/water  cloud  also  has  a  maximum 
concentration  on  the  axis,  and  it  is  likely  that  latent  heat 
released  by  the  condensing  water  vapor  in  this  region  accounts  or 
for  the  secondary  maxima  shown  by  the  particulate  clouds  near  the 
axis  in  the  stratosphere.  When  ambient  moisture  is  included  in 
the  calculation,  however,  the  effects  of  latent  heating  can  be 
much  larger;  Fig.  9  shows  that  for  Case  608,  which  had  a 
moderately  humid  initial  atmosphere  but  was  otherwise  identical 
to  Case  603,  the  plume  exceeds  an  altitude  of  30  km  by 
30  minutes. 
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Figure  9.  Velocity  vector  plot  for  Case  608  at  30  minutes. 
(Includes  ambient  atmospheric  humidity) . 
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(4)  When  the  burning  rate  for  a  10-km  radius  fire  is  changed 
from  the  "fast"  scenario  (Case  502)  to  the  "slow"  scenario 
(Case  504),  peak  temperatures  drop  from  about  1000°K  to  500°K 
(see  Fig.  4),  and  inflow  speeds  in  the  lower  layer  are  reduced  by 
about  a  factor  of  two  (compare  Figs.  5b-c).  Fig.  6c  shows  that 
the  lower  cloud  in  Case  504  stabilizes  at  about  half  the  altitude 
of  the  lower  cloud  in  Case  502  ( 5  km  vs  10  km)  and  spreads 
outwards  at  about  half  the  rate  (reaching  a  radial  distance  of 
30  km  at  2  hours  vs  1  hour) ,  while  the  height  of  the  upper  cloud 
is  reduced  from  25  km  to  15  km.  Similar  results  are  obtained 
from  a  comparison  of  the  multi-phase  Cases  603  ("fast”  burn)  and 
604  ( "slow"  burn) ,  with  the  level  of  the  upper  concentration 
maximum  descending  from  22  km  to  12  km  when  the  burn  rate  is 
reduced,  and  the  lower  maximum  descending  from  12  km  to  5  km 
(compare  Figs.  8  and  10). 


(5)  For  a  fixed  burning  rate,  the  stabilization  height  of  the 
soot  cloud  depends  strongly  on  the  size  of  the  fire.  The  soot 
cloud  for  Case  606,  which  uses  parameters  characteristic  of  the 
largest  World  War  II  firestorms  ("fast”  burn  rate  and  3-km  burn 
radius),  had  a  maximum  altitude  of  14  km  and  was  confined  largely 
to  the  troposphere  (see  Fig.  11).  For  Case  603,  which  used  the 
"fast"  burning  scenario  and  a  radius  of  10  km  (considered  the 
nominal  burn  radius  for  a  1-MT  airburst),  maximum  soot 
concentrations  occur  in  the  stratosphere,  and  the  cloud  top 
extends  to  an  altitude  of  26  km  (see  Fig.  8).  When  the  fire 
radius  is  increased  to  30  km  with  the  same  burn  rate,  the  plume 
exceeds  an  altitude  of  40  km  by  45  minutes  (see  Fig.  12a). 

For  the  inviscid  cases  which  used  the  "fast"  burning 
scenario  and  a  dry  atmosphere,  two  types  of  solution  were 
obtained,  depending  on  the  amount  of  mixing  which  is  allowed  to 
occur  in  the  boundary  layer.  When  the  combustion  heating  is 
confined  to  a  layer  0.1  km  thick,  high  temperatures  result, 
leading  to  overshooting  plumes  and  a  quasi-periodic  solution; 
when  the  combustion  heating  is  "mixed"  over  a  lower  layer  depth 
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Figure  II.  Particulate  clouds  for  Case  606  at  1  hour. 

19 


''  ^  ‘  '  C'  'y^rC'Vv  v  't**f  \*  ,'•  •**•*'*  \«  *.i  ,*  *?  •  *  *  ■-•.-•  v  #l '  -*  *  "  '■■■ 


SOOT  AND  I CE/WATER  CONCENTRATIONS  (mgni/in  ) 


of  1.0  km,  however,  cooler  temperatures  and  less  energetic  plumes 
lead  to  solutions  which  approach  a  steady  state  during  the 
constant  heating  phase. 


Although  in  general  we  do  not  expect  periodic  solutions  from 
non-periodic  heat  sources,  the  fires  which  we  are  simulating  are 
of  unprecedented  size  (with  a  typical  width  more  than  twice  the 
e-folding  depth  of  the  atmosphere),  and  the  oscillatory  behavior 
which  was  obtained  may  have  a  physical  basis.  Mixing  evidently 
plays  a  key  role  in  the  dynamics  of  the  plumes,  however,  and  in 
order  to  resolve  this  question,  a  more  systematic  treatment  of 
its  effects  in  the  model  is  needed. 

Turbulent  Cases 


To  this  end,  a  simplified  turbulence  sub-model  has  been 
developed  and  incorporated  into  the  code.  The  turbulent 
transports  of  heat  and  momentum  are  modeled  according  to 
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where 


aui  3U^ 
ax^  +  axi  ' 


(3) 


9  is  the  potential  temperature,  q  is  the  r.m.s.  turbulent 
velocity,  overbars  denote  ensemble  means,  primes  denote  turbulent 
variations,  and  ilf  f2'  and  Ai  (used  below)  are  turbulent  length 
scales.  Assuming  a  local  balance  between  shear  production, 
buoyant  production,  and  dissipation,  the  above  expressions  may  be 
j  inserted  into  the  turbulent  kinetic  energy  equation  to  give  a 

,  diagnostic  relation  for  qs 
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where  is  the  acceleration  of  gravity  and  9Q  is  a  reference 
temperature . 


We  let  =  hil,  £2  -  A2£,  and  Ki  =  B ,  where  £  is  a 
master  turbulent  length  scale,  and  A2,  A2 ,  and  Bi  are  constants 
known  from  measurements.  To  specify  the  master  turbulent  length 
scale,  we  first  identify  a  natural  length  scale  in  the  flow,  and 
then  determine  the  relationship  between  the  natural  and  turbulent 
length  scales.  In  the  inflow  layer  the  height  above  the  surface 
is  a  natural  length  scale,  and  we  can  use  Blackadar's  [4] 
formula  for  the  turbulent  length  scale 


£  = 
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K  2 

»  KZ  +  ' 

00 


(5) 


where  k  is  von  Karmen's  constant  and  z  is  the  height.  For  large 
values  of  z  the  turbulent  length  scale  asymptotically  approaches 
£w,  which  is  a  prescribed  constant  for  each  case. 

In  order  to  investigate  the  effects  of  parameterized 
turbulence  on  the  solutions.  Cases  702  and  705  were  run  with  the 
same  grid  and  burning  rate  as  Case  502  (see  Table  2),  using  the 
turbulence  sub-model  with  asymptotic  mixing  lengths  of  100  m  and 
25  m,  respectively.  The  main  conclusions  drawn  from  these  two 
cases  were: 


(1)  The  maximum  temperatures  which  developed  (see  Fig.  13)  were 
intermediate  between  those  occurring  in  the  finely  zoned  inviscid 
Case  502  and  the  coarsely  zoned  inviscid  Case  601  (see  Fig.  4). 
Case  705,  which  had  the  smaller  mixing  length  of  the  two 
turbulent  cases,  experienced  less  turbulent  mixing  in  the  inflow 
layer  and  therefore  had  higher  maximum  temperatures  than 
Case  702. 


Table  2.  Turbulent  firestorm  simulation  cases. 


CASE 

BURN 

BURN 

MINIMUM 

SWIRL 

RATE 

RADIUS  (km) 

AZ  (km) 

(m) 

702 

fast 

10 

0.1 

100 

no 

705 

fast 

10 

0.1 

25 

no 

841 

fast 

10 

0.1 

100 

yes 

846 

fast 

10 

0.1 

100—^25 

no 

24 


(2)  The  plume  heights  for  the  two  turbulent  cases  (see 
Figs.  14a-d)  were  also  intermediate  between  those  for  the 
inviscid  Cases  502  (40  km,  not  shown)  and  601  (18  km,  see 
Fig.  7).  The  plume  for  Case  705,  which  had  higher  temperatures 
and  therefore  greater  buoyancy,  reached  a  height  of  over  30  km, 
while  702's  plume  reached  a  maximum  height  of  about  24  km;  the 
soot  tracer  clouds  for  Case  705  were  also  about  3  km  higher  than 
the  corresponding  clouds  for  Case  702  at  45  minutes  (see 
Fig.  15).  Although  both  plumes  showed  considerable  overshoot, 
leading  to  the  formation  of  two  distinct  clouds  of  the  tracers 
released  at  late  times  (note  that  the  tracers  in  Fig.  15  are 
numbered  sequentially  according  to  time-of-release) ,  both  of  the 
turbulent  cases  reached  a  quasi-steady  state  by  30  minutes 
(compare  Figs.  14c-d) . 

These  results  show  that  the  dynamics  of  the  plumes,  and  in 
particular  their  stabilization  heights,  are  sensitive  to  the 
value  of  loo  which  is  chosen.  Based  on  a  formulation  due  to 
Yamada  and  Mellor  [5],  it  can  be  shown  that  i 0 o  should  have  a 
value  between  1/10  and  1/7  of  the  boundary  layer  depth,  and  this 
fact  can  be  used  to  calibrate  the  value  of  the  asymptotic  mixing 
length  to  be  used  for  a  given  problem. 

Approximate  boundary  layer  depths  for  four  of  the  cases  are 
shown  in  the  column  labeled  A  in  the  table  appearing  in  Fig.  16. 
These  boundary  layer  depths  were  estimated  by  taking  the  height  of 
the  25°K  overtemperature  isotherm  at  5  km  radius  at  30  minutes. 

We  picked  this  height  because  this  is  (approximately)  equal  to  the 
cell  depth  in  Case  601.  Since  601  has  no  turbulent  mixing  between 
the  cells,  the  cell  depth  is  the  same  as  the  mixed  or  "turbulent" 
layer  depth  when  the  flow  is  horizontal.  This  is  the  case  for  601 
(see  Fig.  7).  For  Cases  705  and  702,  the  values  of  i»/A  are  1/24 
and  1/8,  respectively.  f«/A  would  have  to  be  smaller  than  1/24 

for  Case  502  and  greater  than  1/8  for  Case  601  to  produce  the 
"observed"  boundary  layer  depths. 


Figure  14a.  Velocity  vector  plots  for  Cases  702  and  705  at  15  minutes. 


Velocity  vector  plots  for  Cases  702  and  705  at  30  minutes 


WHICH  TURBULENCE  FORMULATION  IS  BEST? 


For  boundary  layers  we  can  use 
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where  j  as  0.1 
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Case 

^0 

502 

100  m 

705 

100 

702 

100 
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*  A  is  height  of  +25  K  T1  contour  at  5  km  radius. 
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Figure  16.  Relationship  of  turbulent  mixing  length  to  boundary  layer  depth. 
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Since  Case  702  had  a  value  of  l <»/A  which  lies  between  1/10 
and  1/7,  it  would  appear  to  be  the  most  realistic  of  the 
simulations  discussed  so  far.  The  table  in  Fig.  16  shows  that 
the  boundary  layer  depth  for  this  case  (800  m)  was  closer  to  that 
for  Case  601  (1100  m)  than  to  Case  502  (300  m) ;  the  plume  for 
Case  702,  which  approached  a  steady  state  at  moderate  altitude, 
also  resembles  Case  601 's  plume  (compare  the  left-hand  plots  in 
Figs.  14c-d  with  Figs.  7a-b) .  These  comparisons  show  that  the 
effects  of  parameterized  turbulence  in  a  finely-zoned  grid  are 
similar  to  the  enhanced  mixing  which  results  from  the  use  of  a 
coarse  grid,  and  indicate  that  of  the  inviscid  cases,  the  600- 
series  was  probably  the  more  realistic. 

The  Meteotron  Simulation 


In  order  to  validate  our  model  and  methodology,  the  DICE 
code  was  used  to  simulate  a  600-megawatt  test  fire  conducted  at 
the  Meteotron  research  facility  in  France  on  October  23,  1973 
[6].  The  Meteotron  fire  was  produced  by  the  array  of  oil  burners 
shown  in  Fig.  17,  which  was  simulated  in  our  model  by  a  uniform 
heat  source  of  intensity  150  kW/m2  and  radius  36  m.  The  initial 
atmospheric  temperature  profile  at  the  time  of  the  burn,  shown  in 
Fig.  18,  was  characterized  by  a  sharp  inversion  at  an  altitude  of 
approximately  650  m. 

For  the  large-area  fires  we  have  been  discussing  the  radius 
of  the  burning  area  is  larger  than  the  e-folding  depth  of  the 
atmosphere,  and  the  aspect  ratio  (height  to  width)  of  the 
resulting  plumes  is  of  order  unity.  The  radius  of  the  Meteotron 
fire,  however,  is  much  smaller  than  atmospheric  length  scales, 
and  consequently  we  expect  the  aspect  ratio  of  the  resulting 
plume  to  be  much  greater  than  for  the  firestorm  cases.  In  this 
situation  turbulent  entrainment  into  the  region  of  the  plume 
column  has  a  greater  effect  on  the  buoyancy  than  entrainment  in 
the  inflow  layer.  We  therefore  took  the  natural  length  scale  for 
the  Meteotron  plume  to  be  72  m,  the  width  of  the  burning  region, 
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and  the  asymptotic  length  scale  lM  was  taken  to  be  8  m,  a  value 
which  lies  between  1/10  and  1/7  of  the  natural  length  scale. 

Comparisons  between  the  simulated  and  experimental  Meteotron 
plumes  revealed  the  following  points: 

(1)  During  the  course  of  the  October  23  experiment,  the 
burners  were  run  for  10  minutes,  with  a  "plume  steady  state 
resulting  after  5  minutes;  no  further  increase  of  the  top  plume 
was  visible"  [6].  For  comparison,  the  development  of  the  model 
plume  can  be  most  quickly  visualized  in  Fig.  19,  which  shows  the 
heights  reached  by  tracers  released  at  one  minute  intervals  from 
the  computational  cell  nearest  to  the  center  of  the  fire.  From 
this  figure  it  is  evident  that  the  model  plume  reached  its 
maximum  altitude  between  five  and  six  minutes  into  the  run,  in 
agreement  with  the  experiment. 

( 2 )  The  top  plume  altitude  for  the  experiment  was  reported 
to  be  1260  meters,  while  the  maximum  altitude  of  the  plume  axis 
was  given  as  1050  meters  [6];  the  difference  between  these  two 
heights  is  due  to  the  inclination  of  the  axis  of  the  experimental 
plume  to  the  vertical,  and  may  be  regarded  as  a  rough  measure  of 
"experimental  error."  The  two  reported  heights  have  been  plotted 
with  solid  horizontal  lines  in  Fig.  19,  where  it  can  be  seen  that 
the  model  plume  top  height  was  well  within  the  experimental 
range.  Note  that  both  the  experimental  and  model  plumes 
penetrated  the  inversion  and  rose  a  considerable  distance  into 
the  overlying  stable  layer. 

(3)  Fig.  20  shows  streamlines  of  the  simulated  steady-state 
flow  field  (velocity  is  parallel  to  the  streamlines),  with 
observations  of  the  visible  boundary  of  the  experimental  smoke 
plume  indicated  by  black  circles.  While  turbulent  diffusion 
within  the  plume  causes  the  smoke  to  spread  laterally,  horizontal 
inflow  at  the  edges  of  the  plume  confines  the  smoke  to  the  area 
of  upward  motion.  The  approximate  agreement  between  the  boundary 
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Figure  19. 


Height  reached  by  tracers  released  from  the  computational 
cell  nearest  to  the  center  of  the  fire  for  the  Meteotron 
simulation.  Solid  horizontal  lines  show  experimental  plume 
top  and  plume  axis  heights. 
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of  the  observed  smoke  plume  and  the  edge  of  the  simulated  inflow 
area  indicates  that  the  width  of  plume  has  been  correctly 
predicted . 


(4)  The  averaged  plume  over-temperature  for  eleven 
Meteotron  experiments  is  shown  as  a  function  of  height  by  the 
solid  line  in  Fig.  21,  and  predicted  over-temperatures  for  the 
October  23  simulation  at  one-half  the  observed  plume  radius  are 
indicated  by  black  squares.  In  view  of  the  uncertainties 
involved  in  selecting  a  single  representative  plume  temperature 
at  each  level,  good  agreement  with  the  data  is  indicated. 

These  comparisons  with  the  observed  data  show  that  the  DICE 
simulation  accurately  reproduced  the  dynamics  of  the  October  23 
Meteotron  burn.  Due  to  the  difference  in  aspect  ratios,  physical 
quantities  (e.g.,  mixing  lengths)  cannot  be  scaled  directly  from 
this  case  to  the  firestorm  cases.  In  spite  of  this  difference, 
however,  a  consistent  line  of  physical  reasoning  led  to  the 
selection  of  a  turbulent  length  scale  which  enabled  us  to 
accurately  calculate  the  stabilization  height  of  the  plume  (see 
Fig.  19).  Since  this  experiment  involved  the  penetration  of  a 
thermal  plume  into  an  overlying  stable  layer,  the  success  of  the 
Meteotron  simulation  indicates  that  our  model  is  capable  of 
correctly  simulating  the  penetration  of  a  firestorm  plume  into 
the  stable  layers  of  the  stratosphere. 

The  Effects  of  Swirl 


Some  investigators,  citing  anecdotal  evidence  from  World 
War  II,  have  maintained  that  rotating  flow  is  an  essential 
feature  of  the  firestorm  phenomenon  (e.g.,  [7]).  Others  have 
disputed  this,  and  have  claimed  that  it  is  unlikely  that  large- 
area  fires  could  generate  significant  rotational  velocities  [8]. 
In  order  to  investigate  this  problem,  the  DICE  code  has  been 
extended  to  include  a  tangential  velocity  component;  axial 
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Temperature  (  C) 


Solid  line:  Average  plume  over-temperatures 
for  11  Meteotron  experiments.  Black  squares: 
Model  over-temperatures  calculated  at  one-half 
the  observed  plume  radius. 


symmetry  has  been  preserved,  however,  and  all  quantities  are 
still  functions  of  time,  radius,  and  height  only. 

Case  841  was  run  with  the  extended  version  of  the  code, 
using  the  same  parameters  and  turbulence  sub-model  as  Case  702 
(see  Table  2).  The  initial  swirl  velocity  distribution  for  this 
case  was  given  by 


w  =  2w 


where  x  =  x  /  10  km,  and  w  was  taken  to  be  10  m/sec.  The  profile 
of  the  initial  swirl  velocity  at  the  surface  is  shown  in  Fig.  22; 
in  the  free  atmosphere,  the  initial  swirl  velocity  decreases 
linearly  with  height  from  its  surface  value  to  zero  near  the 
tropopause . 

The  following  conclusions  were  drawn  from  the  swirl  run: 

(1)  By  the  time  the  inflow  layer  is  fully  established  at 

30  minutes,  maximum  swirl  velocities  of  over  200  m/sec  are  found 
in  the  lower  layer  near  the  axis  (see  Fig.  23c,  which  shows  a 
contour  plot  of  the  swirl  speed  at  30  minutes  into  the  run) . 

These  speeds  are  consistent  with  the  conservation  of  angular 
momentum,  since  a  ring  of  air  with  an  initial  tangential  velocity 
of  10  m/sec  at  10  km  radius  would  "spin  up"  to  a  speed  of  about 
300  m/sec  at  a  radial  distance  of  0.33  km  (where  the  innermost 
velocity  gridpoint  is  located),  provided  no  stresses  are  acting. 
While  the  model  has  no  surface  stresses,  turbulent  stresses  from 
the  overlying  air  limit  the  peak  swirl  speeds  in  the  inflow  layer 
to  about  70%  of  the  maximum  possible  value. 

(2)  Theoretical  considerations  lead  us  to  believe  that  "in  a 
nonrotating  convective  cell,  the  vertical  velocity  would 
overshoot  its  equilibrium  altitude ...  in  contrast,  the  swirling 
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Lgure  23c.  Non-swirl  velocity  vectors  and  contours  of  swirl 
speed  for  Case  84 1  at  30  minutes. 


Figure  23d.  Non-swirl  velocity  vectors 


updraft  will  have  an  adverse  pressure  gradient .. .which  not  only 
blocks  the  overshoot,  but  almost  certainly  forces  a  downdraft 
along  the  axis"  ([9],  p.  138).  For  our  baseline  non-rotatinc 
case  (702),  the  plume  shows  considerable  overshoot  by  20  minutes 
(see  Fig.  14b);  for  the  rotating  case  (841),  a  plot  of  the  non¬ 
swirl  velocity  component  at  20  minutes  (see  Fig.  23b)  shows  that 
the  rise  of  the  plume  has  been  substantially  blocked,  and  that  an 
area  of  negative  vertical  velocity  has  indeed  formed  along  the 
axis.  Although  the  swirling  plume  does  penetrate  into  the 
stratosphere  at  later  times  (see  Figs.  23c-d) ,  the  plume  heights 
remain  about  3  km  lower  than  the  corresponding  heights  for  Case 
702  (compare  with  Figs.  14c-d) . 

The  results  of  Case  841  show  that  large  rotational 
velocities  may  develop  in  the  interior  of  a  10-km  radius  fire, 
given  an  initial  ambient  swirl  of  10  m/sec  at  the  fire  edge,  and 
that  the  effect  of  the  swirl  is  to  delay  plume  rise  and  lower  the 
stabilization  height  somewhat.  Interaction  between  the  swirl 
velocity  field  and  the  turbulence  was  not  considered  in  this 
case,  however,  and  it  is  possible  that  the  swirl  may  act  to 
suppress  turbulence,  as  a  stable  density  stratification  does.  As 
was  shown  by  the  comparison  between  Cases  702  and  705  (see 
Figure  14),  this  would  tend  to  intensify  the  plume. 

In  one-dimensional  problems  where  swirl  is  the  only  source 
or  sink  of  turbulent  kinetic  energy  ( TKE ) ,  a  diagnostic  equation 
for  the  TKE  may  be  derived  in  terms  of  the  curvature  Richardson 
number  Rc,  which  is  a  function  of  the  swirl  velocity 
profile  [10].  The  effect  of  curvature  on  the  turbulence  is 
indicated  by  the  stability  function  SC(RC),  shown  in  Fig.  24. 

For  negative  values  of  Sc  the  rate  of  destruction  of  TKE  by  the 
swirl  velocity  field  exceeds  the  rate  of  production,  and  the 
turbulence  is  extinguished. 


In  3-dimensional  cases  where  buoyancy  and  shear  generate 
additional  TKE,  swirl  may  still  exert  a  damping  effect  on  the 
turbulence  or  extinguish  it  completely,  depending  on  the  relative 
magnitudes  of  the  source  and  sink  terms  in  the  TKE  budget 
equation.  Unfortunately  this  equation  is  too  complicated  to  be 
solved  diagnostically  for  the  TKE,  and  to  adequately  treat  the 
general  case,  it  becomes  necessary  to  include  separate  predictive 
equations  in  the  model  for  all  of  the  components  of  the  turbulent 
stress  tensor,  as  well  as  for  the  turbulent  heat  fluxes. 

A  model  developed  by  the  Aeronautical  Research  Associates  of 
Princeton,  Inc.  (ARAP)  currently  has  this  capability  [9].  In  a 
recent  study  of  tornado  dynamics,  two  versions  of  the  ARAP  code 
were  used:  one  employed  a  diagnostic  turbulence  sub-model  similar 
to  that  developed  for  the  DICE  code  (which  omits  the  effect  of 
curvature  on  turbulence),  and  another  used  the  fully  predictive 
sub-model  described  above  (which  includes  the  effects  of 
curvature  on  the  turbulence ) .  Non-dimensional  plots  of  the 
turbulent  eddy  viscosity  field  produced  by  the  diagnostic  and 
predictive  turbulence  models  are  shown  in  Figs.  25a  and  25b, 
respectively.  Near  the  surface  the  turbulence  profiles  are 
similar,  but  at  higher  altitudes  the  predictive  turbulence  model 
produces  much  less  eddy  viscosity,  indicating  that  rotational 
flow  can  have  a  damping  effect  on  turbulence  above  the  boundary 
layer. 

While  the  DICE  code  is  not  yet  capable  of  explicitly 
calculating  the  effects  of  swirl  on  turbulence,  the  likely 
effects  can  be  mimicked  by  artificially  reducing  the  turbulent 
mixing  length  (and  hence  the  turbulence)  in  the  region  above  the 
boundary  layer.  Case  846  was  run  with  the  same  parameters  as 
Case  702  (see  Table  2),  except  that  the  turbulent  mixing  length 
was  reduced  from  100  m  to  25  m  between  the  altitudes  of  0.5  km 
and  1.0  km  (see  Fig.  26).  The  velocity  fields  from  Case  846  at 
20  and  30  minutes  are  shown  in  Figs.  27a-b.  A  comparison  with 
the  corresponding  plots  for  Case  702  (see  Figs.  14b-c)  shows  that 
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decreasing  the  mixing  length  near  the  top  of  the  boundary  layer 
from  its  value  in  Case  702  (100  m)  to  the  value  it  had  in 
Case  705  (25  m)  causes  the  plume  to  rise  by  about  3  km, 
comparable  to  the  decrease  in  plume  height  caused  in  Case  841  by 
swirl  which  does  not  interact  with  turbulence. 

When  the  turbulent  mixing  length  is  changed  from  100  m  to 
25  m  throughout  the  entire  domain,  a  much  greater  increase  in  the 
plume  height  is  obtained,  as  we  saw  in  the  comparison  between 
Cases  702  and  705  (see  Fig.  14).  This  confirms  our  earlier 
observation  that  turbulent  entrainment  within  the  boundary  layer 
is  the  dominant  effect  of  turbulence  in  large-area  fires,  and 
indicates  that  if  swirl  is  to  have  a  significant  effect  on  plume 
stabilization  heights,  it  must  do  so  by  suppressing  turbulence 
within  the  boundary  layer.  Although  present  indications  are  that 
this  effect  would  be  small,  the  definitive  resolution  of  this 
question  awaits  the  modification  of  the  DICE  code  to  predict  the 
full  set  of  turbulent  fluxes,  or  the  application  of  a  model  which 
already  has  that  capability  (such  as  the  ARAP  code)  to  the 
firestorm  problem. 

1.4  CONCLUSIONS  AND  RECOMMENDATIONS. 

The  main  conclusions  which  were  reached  in  the  course  of 
this  study  are  as  follows: 

*  Mixing  plays  a  dominant  role  in  the  dynamics  of  the 

plumes  generated  by  large-area  fires.  When  the  inviscid 
version  of  the  code  is  run  using  the  "fast"  burning 
scenario  and  fine  zoning  (with  a  lower  layer  depth  of 
0.1  km),  high  temperatures  lead  to  the  development  of 
narrow,  overshooting  plumes  and  quasi-periodic  solutions. 
When  the  combustion  heating  is  "mixed"  over  a  lower  layer 
depth  of  1.0  km  in  the  corresponding  coarse-zoned  case, 
cooler  temperatures  lead  to  broader,  less  intense  plumes 
and  a  quasi  steady  state  solution. 
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The  effects  of  mixing  may  be  incorporated  in  a  fine-zoned 
calculation  through  the  use  of  parameterized  turbulence. 
We  used  a  turbulence  sub-model  based  on  the  turbulent 
kinetic  energy  equation,  and  calibrated  the  master 
turbulent  length  scale  through  a  relation  between  the 
turbulent  and  natural  length  scales  proposed  by  Yamada 
and  Mellor  [5].  The  baseline  turbulent  case  (702) 
produced  a  broad,  quasi-steady  state  plume  similar  to 
that  obtained  from  the  coarse-zoned  inviscid  case. 


The  turbulent  version  of  the  model  was  validated  by 
simulating  the  plume  from  an  experimental  fire  conducted 
at  the  Meteotron  research  facility  in  France  on  October 
23,  1973.  Although  the  aspect  ratio  of  the  Meteotron 
plume  was  greater  than  those  of  the  firestorm  plumes,  a 
consistent  line  of  physical  reasoning  led  to  the 
selection  of  a  turbulent  length  scale  which  enabled  us  to 
accurately  calculate  the  stabilization  height  of  the 
plume.  Since  this  experiment  involved  the  penetration  of 
a  thermal  plume  into  an  overlying  stable  layer,  the 
successful  simulation  of  the  Meteotron  plume  indicates 
that  the  DICE  code  is  capable  of  correctly  simulating  the 
penetration  of  a  firestorm  plume  into  the  stable  layers 
of  the  stratosphere. 


The  validation  of  the  turbulent  model  by  the  Meteotron 
simulation,  and  the  similarity  between  the  results  of  the 
baseline  turbulent  case  and  the  single-phase  coarse-zoned 
inviscid  case  (601),  indicate  that  the  soot  and  ice/water 
clouds  generated  in  the  multiphase  coarse-zoned  Cases 
603-606  are  probably  realistic.  These  cases  show  that 
the  stabilization  heights  and  radial  extents  of  the 
particulate  clouds  generated  by  large-area  fires  in  a 
dry,  resting  atmosphere  are  quite  sensitive  to  the  size 
and  burning  rate  of  the  postulated  fire. 
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Latent  heat  released  by  condensing  water  vapor  has  a 
large  effect  on  plume-rise,  since  this  heating  is  usually 
concentrated  at  high  altitudes  near  the  core  of  the 
plume.  When  the  ambient  atmosphere  is  dry,  latent  heat 
released  by  water  vapor  formed  from  the  burning  fuel 
leads  to  the  formation  of  secondary  soot  and  ice/water 
clouds  above  the  central  regions  of  the  primary 
particulate  clouds.  When  the  ambient  atmosphere  is  given 
a  moderate  initial  relative  humidity,  however,  the  effect 
of  latent  heat  release  is  much  more  dramatic,  with  the 
plume  from  a  10-km  "fast"  burning  fire  exceeding  an 
altitude  of  30  Jem  after  30  minutes. 

Swirl  velocities  in  excess  of  200  m/sec  were  obtained  in 
the  interior  of  a  simulated  10-Jan  radius  fire,  given  an 
initial  swirl  velocity  of  10  m/sec  at  the  fire  edge. 

These  large  velocities  are  consistent  with  the 
conservation  of  angular  momentum,  which  holds 
approximately  in  the  boundary  layer  of  the  model  since 
surface  stresses  are  absent  and  stresses  from  the 
overlying  air  are  small.  In  this  simulation,  which  made 
no  allowance  for  the  effects  of  swirl  on  turbulence,  the 
main  effects  of  swirl  were  to  delay  plume  rise  somewhat, 
and  to  lower  the  stabilization  height  of  the  plume  by 
about  3  km . 

When  the  expected  effects  of  swirl  on  turbulence  were 
mimicked  in  a  non-rotating  version  of  the  DICE  code  by 
artificially  reducing  the  turbulent  mixing  length  in  the 
region  above  the  boundary  layer,  the  plume  height  was 
enhanced  by  only  aboc^  3  km,  compared  to  a  control  case 
which  used  a  constant  mixing  length.  When  the  mixing 
length  was  reduced  in  the  boundary  layer  as  well, 
however,  a  much  greater  enhancement  in  plume  rise  was 
obtained.  This  confirms  our  conclusion  that  mixing  in 
the  inflow  layer  is  the  dominant  effect  of  turbulence  in 


large-area  fires,  and  indicates  that  if  swirl  is  to 
dramatically  enhance  the  stabilization  height  of  the 
particulate  clouds  generated  by  large-area  fires,  it  must 
do  so  by  reducing  the  turbulence  in  the  inflow  layer. 

Based  on  these  conclusions,  the  following  recommendations  for 
further  study  are  made: 

•  In  order  to  accurately  assess  the  potential  of  swirl  to 
enhance  plume-rise  from  a  large-area  fire,  it  will  be 
necessary  to  perform  a  firestorm  simulation  which  predicts 
the  full  set  of  turbulent  fluxes,  and  hence  accounts 
correctly  for  the  effects  of  the  swirl  velocity  field  on 
turbulence.  We  therefore  recommend  that  the  swirl  code 
developed  by  the  Aeronautical  Research  Associates  of 
Princeton,  Inc.  (ARAP),  which  already  has  this  capability 
[9],  be  adapted  to  the  firestorm  problem  and  used  to  carry 
out  such  a  study. 

•  The  simulations  presented  in  this  report  showed  that 
persistent,  high-velocity  winds  are  generated  in  the  inflow 
layer  of  a  large-area  fire,  and  it  is  likely  that  these 
winds  will  scour  considerable  quantities  of  debris  from  the 
surface,  particularly  in  areas  of  heavy  blast  damage. 
Scoured  material  which  is  lofted  in  the  plume  may  affect 
the  distribution  of  ice/water  particulates  by  serving  as 
condensation  nuclei,  and  may  affect  the  dynamics  of  the 
plume  through  the  associated  latent  heat  release  as  well  as 
mass  loading  and  drag. 

The  capability  for  predicting  the  input  of  debris  from 
surface  scouring,  as  well  as  its  flow  through  the  grid  in 
discrete  size  groups  which  exchange  heat  and  momentum  with 
the  air,  has  previously  been  implemented  in  the  DICE  code 
in  connection  with  the  study  of  nuclear  burst  phenomena.  A 
new  microphysics  package  has  also  recently  been  implemented 
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in  DICE,  which  predicts  one-parameter  continuous  size 
distributions  for  non-precipitating  cloud  water  and  ice 
hydrometeors  as  well  as  for  rain,  snow,  and  graupel/hail 
particulates  in  each  computational  cell  at  each  time  step. 
If  sub-models  for  the  agglomeration  and  scavenging  of  soot 
are  added  to  these  features,  the  DICE  code  will  be  capable 
of  following  the  evolution  of  the  soot/dust/ice/water  cloud 
from  a  large-area  fire  in  a  fully  interactive  mode. 


We  therefore  recommend  that  this  enhanced  multi-phase 
version  of  the  DICE  code  be  developed,  and  utilized  in 
conjunction  with  a  fine-zoned  grid  and  parameterized 
turbulence,  in  order  to  provide  more  detailed  and  reliable 
estimates  of  the  configuration  and  composition  of  the 
particulate  clouds  generated  by  large-area  fires.  By 
carrying  out  simulations  for  a  variety  of  initial 
temperature  and  humidity  profiles  for  periods  of  up  to 
several  hours  after  the  start  of  the  fire,  considerable 
insight  can  be  gained  into  the  likely  nature  of  post-attack 
environments,  as  well  as  the  potential  for  mass  fires  to 
initiate  the  global  climatic  perturbations  referred  to  as 
"nuclear  winter. " 


SECTION  2 


THE  NUMERICAL  MODEL 


DICE  (an  acronym  for  Dirt  Implicit  Compressible-Fluid 
Eulerian)  is  a  multi-phase,  axisymmetric  numerical  model  based  on 
the  compressible  flow  equations  [3].  The  grid  spacing  can  vary 
both  horizontally  and  vertically,  with  each  cell  representing  an 
average  over  an  annular  region  (see  Fig.  1).  Each  cell  contains 
the  pressure,  density,  and  specific  energy  of  the  air,  which  are 
related  by  the  AFWL  equation  of  state.  Each  cell  also  carries  the 
density  and  specific  energy  for  up  to  nine  additional  groups,  which 
may  be  solid,  liquid,  or  vapor  (e.g.,  soot,  ice/water  particles, 
and  water  vapor) . 


Thermal  interactions  occur  in  a  computational  cell  because 
there  are  (1)  temperature  differences  between  any  of  the  soot- 
water-air  mixture  species,  and  (2)  because  of  material  phase 
changes.  Liquid-vapor  phase  changes  and  the  associated  transfer  of 
latent  heat  are  treated  in  a  time-dependent  fashion.  For  the 
multiphase  firestorm  simulations  discussed  in  this  report  all  of 
the  ice/water  particulates  were  assumed  to  have  a  diameter  of 
10  microns,  and  all  particle  groups  were  flowed  with  the  same 
velocity  as  the  air. 


DICE  utilizes  implicit  finite  differencing,  which  is 
advantageous  in  many  physical  applications  over  the  explicit  finite 
differencing  technique  more  commonly  used  in  wave  propagation 
analyses.  In  the  implicit  method,  the  integration  time  step  is 
inversely  proportional  to  the  maximum  mass  velocities,  whereas  in 
the  explicit  method,  this  step  is  proportional  to  the  reciprocal  of 
the  maximum  wave  velocities.  In  those  problems  where  wave 
velocities  are  much  larger  than  mass  velocities,  an  implicit  method 
integrates  over  a  given  period  of  real  time  in  a  smaller  number  of 
steps  than  is  required  in  an  explicit  method. 
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The  actual  process  of  combustion  is  not  modeled;  instead,  the 
heat  released  by  the  fire  is  added  to  the  lowest  layer  of 
computational  cells  at  a  specified  rate.  The  loss  of  heat  due  to 
thermal  radiation  at  any  location  is  proportional  to  T4  -  Tamb4, 
with  the  constant-  of  proportionality  adjusted  to  make  the  total 
loss  integrated  over  the  domain  equal  to  30%  of  the  rate  of  energy 
input  from  the  fire. 


SECTION  3 


INVISCID  FIRESTORM  SIMULATIONS 


3.1  OVERVIEW. 

This  section  describes  a  series  of  numerical  experiments 
performed  with  an  inviscid  version  of  the  code  (see  Table  1).  The 
results  of  these  cases  illustrate  the  effects  of  variations  in  the 
burning  rate,  burning  area,  computational  zoning,  and  multi-phase 
physics  on  the  predicted  firestorm  environments.  The  initial 
temperature  profile  was  taken  from  the  U.S.  Standard  Atmosphere 
(see  Fig.  2),  for  which  the  tropopause  height  is  11  km.  For  Cases 
301-606,  the  atmosphere  was  assumed  to  be  dry,  while  for  Case  608 
the  initial  relative  humidity  was  assumed  to  be  60%  at  the  surface 
and  to  decrease  linearly  with  altitude  to  zero  at  30  km.  In  Cases 
301-601,  massless  tracer  particles  were  released  in  the  burning 
area  at  regular  intervals;  these  allow  us  to  determine  the  region, 
but  not  the  concentration,  of  the  lofted  soot.  In  Cases  603-608, 

5%  of  the  mass  burned  is  converted  to  soot  and  50%  to  water  vapor 
(see  Table  1),  and  the  soot  and  water  contents  of  the  clouds  are 
calculated  explicitly. 

In  a  study  of  fatalities  in  World  War  II  fires,  Lommasson  and 
Keller  [11]  found  that  these  fires  fell  into  two  broad  classes  (see 
Fig.  28).  Group  fires  had  an  average  burning  rate  of  approximately 
62.5  kw/m2 ,  and  were  characterized  by  relatively  low  fatality  rates 
in  the  population  at  risk.  Firestorms,  which  occurred  less 
frequently,  had  average  burning  rates  of  250  kW/m2,  and  were 
characterized  by  much  higher  fatality  rates  in  the  population  at 
risk.  In  order  to  simulate  these  distinct  types  of  fires,  two 
different  time-varying  burning  rates  were  used  (see  Fig.  3).  In 
the  "fast"  ("slow”)  burning  scenario,  the  burning  rate  rises 
linearly  from  zero  to  a  peak  value  of  250  kW/m2  (62.5  kW/m2 )  in 
15  minutes  (1  hour),  maintains  the  peak  rate  for  30  minutes 


(2  hours),  and  then  decreases  linearly  back  to  zero  in  15  minutes 
(1  hour).  Note  that  for  both  scenarios,  the  total  heat  input  per 
unit  area  integrated  over  the  duration  of  the  fire  is  the  same. 

The  burn  radius  was  taken  to  be  3 ,  10,  or  30  km,  with  the 
combustion  energy  added  uniformly  to  the  lowest  layer  of  cells  over 
the  entire  burning  area.  Cases  301,  502,  and  504  used  a  grid  with 
variable  spacing  and  a  lower  layer  thickness  of  0.1  km,  while  in 
the  600-series  of  cases,  a  uniform  grid  with  1  km  spacing  was  used 
(see  Table  1).  The  radial  and  vertical  boundaries  of  the  grid  were 
placed  at  distances  of  60  km  and  40  km  from  the  origin, 
respectively,  with  continuum  boundary  conditions  allowing  for  the 
flow  of  air  out  of  the  computational  domain  or  the  advection  of  air 
at  ambient  conditions  into  the  domain. 

3.2  FINE-ZONED  CASES. 

In  Case  301,  which  used  the  "fast"  burning  scenario  and  a 
lower  layer  thickness  (where  the  combustion  heating  is  added)  of 
0.1  km,  temperatures  in  excess  of  1200°K  were  generated  (see 
Fig.  29,  which  shows  a  time  history  of  the  temperature  at  three 
locations  within  the  boundary  layer) .  The  large  buoyancy  generated 
by  these  temperatures  led  to  the  formation  of  an  intense,  jet-like 
plume  along  the  axis,  which  developed  vertical  velocities  of  over 
300  m/sec  by  21  minutes  into  the  run  (see  Fig.  30).  These  high 
velocities  caused  the  plume  to  far  overshoot  its  level  of  neutral 
buoyancy,  reaching  altitudes  in  excess  of  40  km,  where  it  became 
overdense  and  fell  back  downwards.  The  downward-moving  jet 
penetrated  near  to  the  surface,  where  it  led  to  the  formation  of  a 
closed  vortex  (see  Fig.  31)  which  eventually  caused  the  upward  jet 
to  weaken  (see  Fig.  32).  The  continued  heating  at  the  surface, 
however,  leads  to  the  formation  of  another  intense,  upward  directed 
jet  along  the  axis  (see  Fig.  33),  and  this  cycle  continues  until 
the  heating  from  the  fire  decreases.  The  resultant  quasi-periodic 
nature  of  the  solution  is  most  evident  in  the  temperature  histories 
(see  Fig.  29),  and  can  also  be  seen  in  the  histories  of  the  inflow 
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Figure  29.  Time  history  of  temperatures  in  the  boundary 
layer  of  Case  301. 
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Figure  30.  Velocity  vector  plot  of  Case  301  at  21  minutes. 
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Figure  33.  Velocity  vector  plot  of  Case  301  at  30  minutes. 
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velocity  at  several  locations  in  the  boundary  layer  (see  Fig.  5a). 
The  oscillatory  behavior  of  the  plume  leads  to  the  formation  of  two 
distinct  regions  of  lofted  soot,  as  indicated  by  the  positions 
reached  at  one  hour  by  massless  tracers  released  every  five  minutes 
from  within  the  burning  area  (see  Fig.  6a).  One  region  near  the 
axis  at  an  altitude  of  about  30  km  is  associated  with  the 
overshooting  phase  of  the  plume,  and  the  other,  extending  radially 
outwards  to  a  distance  of  30  km  near  the  altitude  of  the 
tropopause,  shows  the  location  of  the  outflow  area. 

Case  502  used  the  same  parameters  as  Case  301,  except  that  the 
grid  spacing  in  the  radial  direction  was  twice  as  large.  In  spite 
of  this  difference,  the  flow  field  for  Case  502  exhibits  upward  and 
downward  jets  and  closed  vortices  similar  to  those  obtained  in  Case 
301,  as  shown  by  velocity  vector  plots  of  the  two  cases  at  25 
minutes  (see  Figs.  34a-b) .  Like  Case  301,  Case  502  exhibited  a 
quasi-periodicity  of  about  ten  minutes,  as  can  be  most  readily  seen 
in  the  maximum  temperature  history  for  that  case  (see  Fig.  4). 

These  temperatures  were  slightly  less  than  the  maximum  temperatures 
for  Case  301  (compare  with  Fig.  29),  since  Case  301  had  finer 
radial  zoning  near  the  axis  and  was  therefore  better  able  to 
resolve  the  hottest  regions  of  the  plume.  Consequently  the  soot 
tracer  clouds  for  Case  502,  while  showing  two  regions  of  maximum 
concentration  as  they  did  for  Case  301,  stabilized  at  slightly 
lower  altitudes  (compare  Figs.  6a-b) . 

Case  504  was  identical  to  Case  502,  except  that  the  "slow" 
burning  scenario  was  used.  The  factor  of  4  reduction  in  the 
burning  rate  caused  peak  inward  velocities  at  various  ranges  near 
the  ground  to  be  reduced  by  about  a  factor  of  2.  Figs.  5b-c 
compare  radial  velocity  histories  near  the  ground  at  ranges  between 
8.5  and  20  km.  At  the  8.5  km  range,  peak  inward  winds  of  about 
40  m/sec  (90  mph)  are  reduced  to  20  m/sec  for  the  high  and  low 
burning  rate  cases,  respectively.  Time  histories  of  the  maximum 
inward  wind  speeds,  which  are  reached  near  the  axis  of  the  fire, 
are  shown  for  Cases  502  and  504  in  Fig.  5d.  Peak  velocities  retch 


Velocity  vector  plot  of  Case  301  Figure  34b.  Velocity  vector  plot  of  Case  302 

at  25  minutes.  at  25  minutes. 


100  m/sec  (220  mph)  in  the  firestorm  case,  versus  50  m/sec 
(110  mph)  for  the  group  fire  simulation.  Buoyancy  forces  are  the 
primary  cause  for  the  large  wind  speeds.  One  measure  of  this  force 
is  the  maximum  temperature  in  the  burning  region.  The  maximum 
temperature  histories  versus  time  for  these  two  cases  may  be  seen 
in  Fig.  4.  Maximum  temperatures  reach  about  1000°K  compared  to 
500°K  in  Cases  502  and  504,  respectively. 

The  factor  of  4  reduction  in  burning  rate  causes  the  maximum 
stabilization  altitudes  for  lofted  soot  particles  to  be  reduced  by 
almost  a  factor  of  2.  Figs.  6b-c  show  a  comparison  of  the 
positions  of  lofted  soot  tracers  when  the  cloud  has  expanded  about 
30  km  in  radius;  this  radial  expansion  takes  two  hours  in  the 
reduced  burning  rate  case  compared  to  one  hour  in  the  higher  rate 
case.  There  are  two  distinct  regions  of  lofted  particles  in  both 
cases.  The  higher  region,  near  the  axis,  stabilized  at  a  maximum 
altitude  of  about  25  km  (Case  502)  compared  to  15  km  (Case  504). 

The  lower  region  of  lofted  particles  forms  a  radially  expanding 
volume  with  a  mean  altitude  of  about  10  km  versus  5  km  in  the  high 
versus  low  burning  rate  cases,  respectively. 

3.3  COARSE-ZONED  CASES. 

The  remaining  cases  to  be  discussed  in  this  section  (the  600- 
series  in  Table  1)  used  a  uniform  grid  with  cell  dimensions  equal 
to  1  km.  Since  the  combustion  heating  is  added  in  the  lowest  layer 
of  cells,  the  energy  input  from  the  fire  in  these  cases  is 
effectively  "mixed"  over  the  lowest  1  km  of  the  computational 
domain,  thereby  crudely  simulating  the  effects  of  turbulence  in  the 
inflow  layer. 

Case  601  is  identical  to  Case  502,  except  for  the  change  in 
zoning.  Maximum  temperatures  for  this  case  are  much  cooler  than 
for  502,  and  reach  a  steady  state  by  30  minutes  (see  Fig.  4).  The 
plume  for  case  601  also  reaches  a  steady  state  by  30  minutes 
(compare  Fig3.  7a-b),  and  exhibits  neither  the  intense,  narrow  jets 


nor  the  closed  vortices  shown  by  Cases  301  and  502.  The  maximum 
altitude  reached  by  601 's  plume  (18  km)  is  also  lower  than  the 
maximum  reached  by  Case  502  (40  km,  not  shown).  Since  the  plume 
for  Case  601  shows  little  tendency  to  overshoot,  the  associated 
tracer  cloud  (Fig.  6d)  is  concentrated  near  the  tropopause,  and 
does  not  show  a  second  maximum  at  higher  levels  as  in  Cases  301 
and  502,  which  had  finer  zoning  and  therefore  hotter  driving 
temperatures . 

In  Cases  603-608,  5%  of  the  mass  burned  is  converted  to  soot, 
50%  is  converted  to  water  vapor,  and  concentrations  of  soot,  water 
vapor,  liquid  water,  and  ice  are  predicted  by  the  model.  Although 
phase  changes  of  water  are  included,  precipitation  is  not;  most  of 
the  water  vapor  condenses  into  liquid,  and  eventually  freezes  into 
suspended  ice  particles  as  the  cloud  rises  and  cools  adiabatically . 

Case  603  differs  from  Case  601  only  in  the  inclusion  of  soot 
and  water  species  in  the  calculation.  The  soot  cloud  predicted  for 
Case  603  at  one  hour  (see  Fig.  8)  shows  a  maximum  concentration  in 
the  lower  stratosphere,  as  does  the  tracer  cloud  for  Case  601 
(see  Fig.  6d) .  The  configuration  of  the  ice/water  cloud  predicted 
at  one  hour  for  Case  603  resembles  that  of  the  soot  cloud  (see 
Fig.  8),  with  a  "shield"  of  particulates  extending  outwards  above 
the  tropopause  to  a  radial  distance  of  about  40  km  after  one  hour. 
Concentration  levels  for  the  ice/water  cloud  are  about  an  order  of 
magnitude  greater  than  those  for  the  soot  cloud,  since  ten  times  as 
much  fuel  mass  is  converted  into  water  vapor  as  is  converted  into 
soot  (see  Table  1). 

Both  of  the  particulate  clouds  for  Case  603  show  considerable 
vertical  dispersion,  compared  to  the  relatively  thin  outflow  layer 
shown  by  the  tracer  cloud  for  Case  601  (see  Fig.  6d);  the  ice/water 
cloud  is  confined  to  altitudes  above  6  km,  however,  since  water 
vapor  below  this  altitude  has  not  yet  undergone  sufficient  cooling 
to  condense.  Both  clouds  show  maximum  concentration  near  the  axis, 
indicating  that  most  of  the  condensation  of  water  vapor  occurs  in 
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this  region.  Although  the  latent  heat  released  by  the  condensing 
water  vapor  amounts  to  only  5%  of  the  total  combustion  energy 
released  by  the  fire,  the  concentration  of  this  heating  in  a 
relatively  small  area  in  the  core  of  the  plume  causes  it  to 
overshoot  slightly,  compared  to  the  plume  for  Case  601;  this  in 
turn  leads  to  the  formation  of  secondary  maxima  in  both  of  the 
particulate  clouds  at  altitudes  above  20  km  near  the  axis. 

Case  604  is  the  same  as  Case  603  except  that  the  "slow" 
burning  scenario  is  used.  Fig.  10  shows  that  in  this  case  also, 
the  release  of  latent  heat  by  condensing  water  vapor  near  the  axis 
leads  to  a  double-lobed  soot  cloud  (note  that  the  lower  ice  cloud 
is  missing,  since  water  vapor  which  detrains  below  5  km  has  not  yet 
undergone  sufficient  cooling  to  condense).  As  for  Case  504  the 
slower  burning  rate  leads  to  lower  plume  temperatures,  and  the 
reduced  buoyancy  which  results  causes  the  plume  to  stabilize  at  a 
lower  height.  The  upper  concentration  maximum  descends  from  22  km 
in  Case  603  to  12  km  in  Case  604,  and  the  lower  concentration 
maximum  descends  from  12  km  to  5  km,  as  the  burning  rate  is  reduced 
by  a  factor  of  4 . 

Case  605  is  also  identical  to  Case  603,  except  that  the 
burning  radius  for  605  was  increased  from  10  to  30  km;  thus 
although  both  cases  have  the  same  energy  release  rate  per  unit 
area,  the  total  energy  release  for  Case  605  is  almost  an  order  of 
magnitude  greater  than  603 's.  For  the  first  30  minutes  the  main 
body  of  605 's  plume  does  not  penetrate  above  18  km  (see  Fig.  12a), 
and  the  soot  and  ice/water  clouds  at  that  time,  while  containing 
more  total  mass  than  the  clouds  for  Case  603  (since  more  fuel  has 
been  burned),  are  also  confined  below  that  altitude  (see  Fig.  12b). 
By  45  minutes  into  the  run,  however,  as  the  inflow  from  the  larger 
heating  area  reaches  the  axis,  an  intense  plume  penetrates  above 
40  km  altitude  (see  Fig.  12a). 

Case  606  is  the  same  as  Case  603  except  that  the  burn  radius 
is  reduced  from  10  km  to  3  km.  In  this  case,  as  in  the  "slow" 


burning  scenario  (Case  604),  most  of  the  soot  is  confined  to  the 
troposphere,  although  some  portions  of  the  soot  cloud  penetrate 
into  the  stratosphere  near  the  axis  (see  Fig.  11).  The  ice/water 
cloud  is  also  largely  confined  to  the  troposphere,  although  the 
altitude  of  maximum  concentration  is  somewhat  higher  than  that  of 
the  soot  cloud,  since  water  vapor  in  the  lower  troposphere  has  not 
been  cooled  enough  to  condense.  The  lateral  spreading  rate  of  the 
particulate  clouds  are  also  greatly  reduced;  the  ice/water  cloud, 
in  particular,  reaches  a  radial  extent  of  only  6  km  after  one  hour, 
compared  with  30  km  in  Case  603,  which  had  a  burn  radius  of  10  km 
( see  Fig.  8 ) . 

Case  608  is  the  same  as  Case  603  except  that  the  ambient 
atmosphere  is  moist,  with  the  initial  relative  humidity  varying 
linearly  from  60%  at  the  surface  to  zero  at  an  altitude  of  30  km. 

In  Case  603,  latent  heat  released  by  water  vapor  which  formed 
during  combustion  of  the  fuel  led  to  a  moderate  increase  in  plume 
height  (26  km  compared  to  18  km  for  Case  601).  When  ambient 
moisture  from  a  moderately  humid  atmosphere  is  added,  however,  the 
effects  of  latent  heat  release  are  much  greater,  with  the  plume  for 
Case  608  passing  an  altitude  of  30  km  by  30  minutes  into  the  run 
( see  Fig.  9 ) . 

3.4  DISCUSSION. 

Cases  606,  603,  and  605,  all  of  which  used  the  "fast"  burning 
scenario,  illustrate  the  effects  of  increasing  the  burn  radius  from 
3  to  10  to  30  km.  In  Case  606,  for  which  the  burning  area  was 
similar  to  that  which  characterized  the  largest  firestorms 
occurring  in  the  Second  World  War,  most  of  the  soot  is  confined  to 
the  troposphere,  with  the  cloud  top  reaching  14  km.  When  the 
burning  radius  is  increased  to  10  km  (considered  the  nominal  burn 
radius  for  a  1-MT  airburst),  the  height  of  the  soot  cloud  is  raised 
considerably,  with  the  maximum  concentration  now  occurring  above 
the  tropopause  (see  Fig.  8),  and  the  cloud  top  reaching  26  km. 


Increasing  the  burn  radius  to  30  km  results  in  another  large 
increase  in  the  plume  height  (see  Fig.  12a). 

The  "fast"  burning  scenario  is  intended  to  simulate  the 
burning  rate  from  World  War  II  firestorms,  which  occurred  in  old, 
high-density  city  centers.  Use  of  the  "slow"  burning  scenario, 
which  may  be  more  representative  of  modern,  low-density  cities, 
results  in  a  considerable  lowering  of  the  estimated  cloud  height 
for  a  10  km  radius  fire  (compare  Case  603  with  604,  and  502  with 
504).  These  cases  show  that  predicted  cloud  stabilization  heights, 
and  in  particular  the  degree  of  penetration  of  soot  into  the 
stratosphere,  are  dependent  on  the  radius  and  burn  rate  of  the 
postulated  fire.  Case  608  also  showed  that  the  plume  height  may  be 
quite  sensitive  to  the  humidity  of  the  ambient  atmosphere  at  the 
time  of  the  fire. 

For  the  cases  which  used  the  "fast"  burning  scenario  and  a  dry 
atmosphere,  two  types  of  solution  were  obtained,  depending  on  the 
amount  of  mixing  which  is  allowed  to  occur  in  the  boundary  layer. 
When  the  combustion  heating  is  confined  to  a  layer  0.1  km  thick, 
high  temperatures  result,  leading  to  overshooting  plumes  and  a 
quasi-periodic  solution;  when  the  combustion  heating  is  "mixed” 
over  a  lower  layer  depth  of  1.0  3cm,  however,  cooler  temperatures 
and  less  energetic  plumes  lead  to  solutions  which  approached  a 
steady-state  during  the  constant  heating  phase.  Although  in 
general  we  do  not  expect  periodic  solutions  from  non-periodic  heat 
sources,  the  fires  which  we  are  simulating  are  of  unprecedented 
size  (with  a  typical  width  more  than  twice  the  e-folding  depth  of 
the  atmosphere),  and  the  oscillatory  behavior  which  we  obtained  may 
have  a  physical  basis.  Mixing  evidently  plays  a  key  role  in  the 
dynamics  of  the  plumes,  however,  and  in  order  to  resolve  this 
question,  a  more  systematic  treatment  of  its  effects  in  the  model 
is  needed.  To  this  end,  a  simplified  turbulence  model  has  been 
developed  and  incorporated  into  the  code. 


SECTION  4 
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THE  EFFECTS  OF  PARAMETERIZED  TURBULENCE 


4.1  FORMULATION  OF  THE  MODEL. 


Turbulent  motions,  which  involve  a  wide  range  of  length  scales 
and  are  inherently  three-dimensional,  are  not  calculated  explicitly 
in  our  model.  Instead,  the  turbulence  is  “averaged  out"  by 
applying  an  averaging  operation  to  the  original  governing 
equations.  For  example,  the  governing  mass  continuity  and  momentum 
equations  for  a  perfect  incompressible  fluid  are 


where  Uj  is  the  j-th  component  of  the  velocity  vector,  p  is 
pressure,  p  the  density,  and  Fj  the  net  body  force  acting  on  the 
fluid.  Divide  the  variables  into  a  mean  part  (denoted  by  an 
over bar)  and  a  deviation  from  the  mean: 
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P  =  P  +  P' 
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Substitute  these  into  Eqn.  (1)  and  average  to  get  the  mean 
governing  equations: 
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Note  the  additional  term  in  Eqn.  (9).  This  represents  the 
turbulent  momentum  transport,  and  adds  three  unknowns  to  the  four 
we  already  have.  Since  we  still  have  only  four  equations  (one 
continuity  and  three  momentum),  the  set  is  not  closed.  To  close 
the  set,  we  must  somehow  specify  the  turbulent  fluxes.  This  is 
the  turbulence  closure  problem. 

Many  turbulence  closure  models  have  been  developed.  Some 
assume  that  turbulent  transport  is  like  molecular  diffusion. 
Others  use  equations  such  as  the  turbulent  kinetic  energy 
equation  in  their  models,  and  take  the  buoyant  generation  of 
turbulent  kinetic  energy  into  account.  How  do  we  choose  a 
turbulence  model  that  is  appropriate  for  the  flow  of  interest? 

Our  philosophy  is  to  use  the  simplest  model  that  includes  the 
processes  that  create  or  destroy  turbulent  kinetic  energy  in  the 
flow  of  interest. 

We  will  use  the  Boussinesq  equations  in  Cartesian 
coordinates.  The  averaged  continuity,  momentum,  and 
thermodynamic  equations  are 
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Here  9  =  T(p0/p)0*277  is  the  potential  temperature,  where  T  is 
the  temperature.  PQ  is  a  constant  reference  pressure,  60  is 
a  constant  reference  temperature,  and  U ^'9'  is  the  turbulent 
flux  of  potential  temperature.  We  will  specify  the  unknown 
turbulent  covariances  in  terms  of  the  mean  fields  as  follows: 
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where  q  is  the  r.m.s.  turbulent  velocity  (a  measure  of  the 
turbulent  kinetic  energy) ,  and  and  l2  are  turbulent  length 
scales.  These  models  are  strictly  valid  only  in  the  limit  of 
isotropic,  nearly  homogeneous  turbulence. 


To  close  these  equations  (treating  Ui  and  6  as  known)  we 
need  an  equation  for  q2 .  This  is  the  turbulent  kinetic  energy 
equation: 


tKuruc 

i  i  j 


2 


'au±  au^ 

axT  +  ax. 

j 


+  2 

eo 


uisi 


2e  (12) 


The  first  term  on  the  right  represents  the  turbulent  transport, 
the  second  the  shear  production,  the  third  the  buoyant 
production,  and  the  last  is  molecular  dissipation.  We  assume 
that  production  and  dissipation  balance: 


0  =  S.P.  +  B.P.  -  D, 


(13) 


use  the  above  models  for  the  fluxes,  and  model  dissipation  by 


e  =  q  /A. . 


(14) 


The  result  is 
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where 
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Let  i  1  =  Axi  ,  i2  =  ^2  >  an<^  Ai  =  B^  ,  where  i  is  a  master 
turbulent  length  scale,  and  Alf  A2 ,  and  Bx  are  constants  known 
from  measurements. 

We  still  must  specify  the  master  turbulent  length  scale.  To 
do  this,  we  first  identify  a  natural  length  scale  in  the  flow. 
Then  we  must  determine  the  relationship  between  the  natural  and 
turbulent  length  scales.  If  we  use  measurements  from  actual 
turbulent  flows  to  determine  this  relationship,  we  must  then 
decide  whether  this  relationship  is  likely  to  hold  in  the  flow  we 
wish  to  simulate. 

A  plume  has  three  regions:  the  inflow  layer,  the  plume 
column,  and  the  outflow  layer.  Each  has  a  natural  length  scale. 
The  inflow  layer's  depth  is  a  natural  length  scale,  as  is  the 
width  of  the  plume  column  arid  the  depth  of  the  outflow  layer. 

What  are  the  relationships  between  these  natural  length  scales 
and  the  turbulent  length  scales?  The  inflow  layer  is  a  boundary 
layer,  so  we  can  use  Blackadar's  [4]  formula: 
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where  k  is  von  Karmen's  constant,  z  is  the  height  above  the 
surface,  and  £«  is  related  to  the  depth  of  the  boundary  layer. 
Yamada  and  Mellor  [5]  related  i*  directly  to  the  turbulence 
profile: 

f* 

i„=  2a  d*  (18) 
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where  a  =  0.1.  We  could  use  the  turbulence  model-derived  values 
of  q,  the  r.m.s.  turbulent  velocity,  to  estimate  l».  This  could 
be  done  at  each  radial  location  every  time-step  while  the 
simulation  proceeds.  We  didn't  do  this.  Instead,  we  simply 
specified  for  each  case. 

Far  from  a  boundary,  the  turbulent  length  scale  is  related 
directly  to  the  width  of  the  turbulent  region,  and  we  could  use 
the  formula  above  to  relate  i<*>  to  the  turbulence  profile  across 
the  plume  column  or  outflow  layer  at  each  time  step.  Since  the 
aspect  ratio  of  the  firestorm  plume  is  only  of  order  unity, 
however  (see  Fig.  7),  it  will  have  a  large  volume  to  surface  area 
ratio.  We  expect  the  effects  of  turbulent  entrainment  to  be 
small  in  the  region  of  the  plume  column,  therefore,  and  so  for 
simplicity,  Eqn.  (17)  was  used  to  specify  the  turbulent  mixing 
length  once  and  for  all  throughout  the  computational  domain. 

4.2  CASES  702  AND  705. 

These  cases  were  similar  to  Case  502,  in  that  they  used  the 
"fast"  burning  scenario  and  the  same  grid;  the  turbulence  model 
was  included,  however,  with  the  master  turbulent  length  scale, 
loo,  set  equal  to  100  meters  for  Case  702  and  25  meters  for 
Case  705  (see  Table  2). 

Fig.  14  shows  the  time  evolution  of  the  velocity  fields.  At 
15  minutes,  the  inflow  layers  are  nearly  fully  developed  and 
contain  peak  winds  of  more  than  50  m/s,  but  the  plumes  are  just 
starting  to  develop.  During  the  next  five  minutes,  the  plumes 
shoot  upwards  at  50  m/s  (Case  702)  and  80  m/s  (Case  705)  and 
reach  21  km  (Case  702)  and  30  km  (Case  705).  By  30  minutes,  the 
plumes  are  nearly  steady,  and  reach  maximum  heights  of  24  km 
(Case  702)  and  more  than  30  km  (Case  705). 
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After  ascending  in  the  plumes,  the  air  descends  and 
oscillates  about  its  non-buoyancy  level  slightly  before 
stabilizing  in  the  outflow  layer.  This  layer  extends  from  about 
12  to  15  km  altitude  in  Case  702  =  100  m)  and  from  about  13 

to  18  km  altitude  in  Case  705  (2»  =  25  m) ,  as  shown  in 
Fig.  14d.  We  also  released  tracer  particles  at  the  surface 
within  the  burning  region  at  five  minute  intervals.  Their 
locations  should  indicate  the  top  of  the  stabilization  or  outflow 
layer  of  the  plume,  since  they  will  follow  the  surface  air,  which 
(we  will  see)  is  the  warmest  and  therefore  most  buoyant.  The 
tracer  particle  locations  at  45  minutes  are  shown  in  Fig.  15; 
they  are  labeled  sequentially  according  to  their  time-of-release . 
The  tracer  particles  indicate  that  the  outflow  layer  top  is  about 
15  km  in  702  and  from  15  to  19  km  in  Case  705. 

Fig.  35(a)  compares  the  temperature  departure  from  the 
initial  conditions  at  30  minutes  for  Cases  702  and  705;  Fig. 

35(b)  is  the  same  but  shows  only  the  lowest  and  innermost  10  km. 
Since  the  temperature  departure  is  proportional  to  the  buoyancy 
force,  differences  in  the  temperature  fields  are  related  to  the 
differences  in  the  plume  dynamics.  The  plume  in  Case  705 
(i»  =  25  m)  reaches  higher  than  does  Case  702 's  plume 
(i»  =  100  m) .  Fig.  35(a)  shows  that  at  all  altitudes  below 
21  km,  Case  705 's  plume  contains  warmer  air  at  the  axis  than  does 
Case  702 's  plume,  which  is  consistent  with  its  greater  height. 
Figs.  35(b)  and  13  show  that  the  inflow  layer  became  about  100°K 
warmer  in  Case  705  than  in  Case  702.  This  difference  is  due  to 
greater  turbulent  mixing  in  Case  702 's  inflow  layer. 

One  might  ask  whether  the  inflow  is  more  rapid  in  Case  702 
than  in  Case  705,  which  would  allow  less  time  for  the  air  to  heat 
up.  Fig.  36  shows  that  this  isn't  the  case;  the  maximum  inflow 
speed  is  about  50  m/s  in  Case  702  and  60  m/s  in  Case  705. 

For  comparison,  temperature  departure  fields  at  30  minutes 
for  the  inviscid  Cases  502  ( A Z 0  =  100  m)  and  601  ( A Z 0  =  1000  m) 
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Figure  36.  Maximum  inflow  velocity  for  Cases  702  and  705 
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are  shown  in  Fig.  37.  The  most  striking  difference  between  the 
temperature  departure  fields  for  these  two  cases  may  be  seen  in 
the  inflow  layer,  where  over-temperatures  reach  400°K  in 
Case  502,  but  only  50°K  in  Case  601.  The  temperature  departure 
fields  for  the  turbulent  Cases  702  and  705  (see  Fig.  35b)  are 
intermediate  between  those  for  601  and  502.  This  means  that 
larger  grid  cells  in  the  lowest  layer  increase  the  mixing  in  the 
inflow  layer.  This  occurs  because  distributing  the  burning 
through  the  lowest  layer  of  large  cells  is  equivalent  to 
distributing  the  burning  through  several  thinner  cell  layers. 

From  a  comparison  of  Cases  702  and  705,  it  is  evident  that 
the  plume  dynamics,  and  in  particular  the  stabilization  height, 
depend  on  the  asymptotic  turbulent  length  scale  (i»)  which  is 
specified.  How  do  we  decide  which  length  scale  is  best? 

Fig.  16  shows  how  £«,  as  defined  in  Eqn.  18,  depends  on  the 
depth  (A)  of  the  boundary  layer  for  three  different  profiles  of 
q.  For  a  uniform  distribution  of  q  with  height,  i<»  is  a  tenth  of 
the  depth  of  the  boundary  layer,  while  for  a  Gaussian 
distribution  of  q,  Jt„  is  0.14  times  the  half-width  of  the 
profile.  Thus,  we  expect  that  the  most  realistic  simulation  will 
be  the  one  in  which  the  prescribed  loo  is  1/10  to  1/7  of  the 
turbulent  boundary  layer  depth. 

We  estimated  the  turbulent  boundary  layer  depth  to  be  the 
height  of  the  25°K  temperature  departure  isotherm  at  5  km  radius 
at  30  minutes .  We  picked  this  height  because  this  is 
(approximately)  equal  to  the  cell  depth  in  Case  601.  Since  601 
has  no  turbulent  mixing  between  cells,  the  cell  depth  is  the  same 
as  the  mixed  or  "turbulent"  layer  depth  when  the  flow  is  purely 
horizontal  between  5  and  10  km.  This  is  the  case  for  601;  see 
Fig.  7. 


The  estimated  turbulent  boundary  layer  depths  are  shown  in 
Fig.  16  under  the  column  heading  A.  They  range  from  300  m  in 
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Figure  37.  Temperature  departure  at  30  minutes  for  Cases  502  and  601. 


Case  502  to  1100  m  in  Case  601.  Also  listed  are  £«/ A.  This  is 
1/8  for  Case  702  and  1/24  for  705.  £»/A  would  be  have  to  be 
smaller  than  1/24  for  Case  502  and  larger  than  1/8  for  Case  601 
to  produce  the  "observed"  boundary  layer  depths. 

4.3  DISCUSSION. 

Since  Case  702  had  a  value  of  £«/A  which  lies  between  1/10 
and  1/7,  it  would  appear  to  be  the  most  realistic  of  the 
simulations  discussed  so  far.  The  table  in  Fig.  16  shows  that 
the  boundary  layer  depth  for  this  Case  (800  m)  was  closer  to  that 
for  Case  601  (1100  m)  than  to  Case  502  (300  m) ,  and  the  plume  for 
Case  702,  which  approached  a  steady-state  at  moderate  altitude, 
also  resembles  Case  601's  plume  (compare  Figs.  14cd  and  7ab) . 

This  shows  that  the  effects  of  parameterized  turbulence  in  a 
finely- zoned  grid  are  similar  to  the  enhanced  mixing  which 
results  from  the  use  of  a  coarse  grid,  and  indicates  that  of  the 
inviscid  cases,  the  600-series  was  probably  the  more  realistic. 

The  results  of  Case  702  indicate  that  most  of  the 
contaminants  from  a  10-km  radius  fire,  with  a  peak  burning  rate 
of  250  kw/m2,  would  be  expected  to  detrain  at  levels  well  within 
the  stratosphere  (see  Fig.  15).  How  much  credibility  can  be 
attached  to  results  such  as  these?  In  order  to  be  confident  that 
the  dynamics  of  the  thermal  plumes  are  being  accurately 
represented,  the  numerical  models  used  should  be  verified  with 
experimental  data.  Experiments  on  the  firestorm  scale  are 
impractical,  of  course,  but  verification  with  smaller-scale  test 
fires  may  partially  satisfy  this  requirement.  One  such  study 
using  the  DICE  code  is  discussed  in  the  following  section. 


SECTION  5 


THE  METEOTRON  EXPERIMENT 


5.1  OVERVIEW. 


We  used  our  model  to  simulate  a  600-Megawatt  test  fire 
conducted  at  the  Meteotron  facility  of  the  Centre  de  Recherches 
Atmospheriques ,  France,  on  October  23,  1973.  The  experimental 
plume  at  the  Meteotron  facility  was  produced  by  the  array  of  oil 
burners  depicted  in  Fig.  17.  Our  model  did  not  attempt  to 
resolve  individual  flames,  but  assumed  a  steady,  uniform  heat 
input  of  150  kW/m2  into  the  lowest  layer  of  computational  cells, 
over  a  circular  area  of  radius  36  m.  The  computational  domain 
measured  1  km  radially  by  2  km  vertically,  with  continuum 
boundary  conditions  to  allow  for  the  flow  of  air  out  of  or  into 
the  domain. 

5.2  INITIAL  CONDITIONS. 

The  ambient  atmospheric  temperature  profile  measured  at  the 
time  of  the  burn  is  shown  in  Fig.  18.  The  lapse  rate  between  the 
surface  and  650  m  was  about  equal  to  the  dry  adiabatic  lapse  rate 
of  -1°C/100  m,  indicating  a  layer  of  neutral  stability.  At  650  m 
there  was  a  sharp  inversion,  with  temperatures  continuing  to 
increase  up  to  a  height  of  1250  m,  indicating  a  layer  of  high 
static  stability. 

A  contour  plot  of  the  initial  potential  temperature  field  is 
shown  in  Fig.  38a.  In  this  representation,  the  layer  of  neutral 
stability  appears  as  a  region  of  constant  potential  temperature; 
the  overlying  stable  layer,  characterized  by  a  rapid  increase  of 
potential  temperature  with  height,  appears  as  a  region  of  closely 
packed  isentropes  (lines  of  constant  potential  temperature).  An 
isentropic  plot  of  the  U.S.  Standard  Atmosphere  is  shown  for 


Figure  38a.  Initial  potential  temperature  Figure  38b.  Potential  temperature  field  for 

field  for  the  Meteotron  the  U.S.  standard  atmosphere. 

experiment. 


comparison  in  Fig.  38b  (note  the  different  scaling).  In  this 
figure  the  troposphere  appears  as  a  region  of  low  static 
stability  (weak  potential  temperature  gradient),  while  the 
stratosphere  appears  as  a  region  of  high  static  stability  (strong 
potential  temperature  gradient).  These  figures  show  that  the 
structure  of  the  lower  atmosphere  for  the  October  23  Meteotron 
burn  resembled  a  scaled-down  version  of  the  troposphere- 
stratosphere  system. 

For  the  large-area  fires  we  have  been  discussing  the  radius 
of  the  burning  area  is  greater  than  the  e-folding  depth  of  the 
atmosphere,  and  the  aspect  ratio  (height  to  width)  of  the 
resulting  plumes  is  of  order  unity.  The  radius  of  the  Meteotron 
fire,  however,  is  much  smaller  than  atmospheric  length  scales, 
and  consequently  we  expect  the  aspect  ratio  of  the  resulting 
plume  to  be  much  larger  than  for  the  firestorm  cases.  In  this 
situation  turbulent  entrainment  into  the  region  of  the  plume 
column  has  a  greater  effect  on  the  buoyancy  than  entrainment  in 
the  inflow  layer.  We,  therefore,  took  the  natural  length  scale 
for  the  Meteotron  plume  to  be  72  m,  the  width  of  the  burning 
region,  and  the  turbulent  length  scale  SL »  was  taken  to  be  8  m,  a 
value  which  lies  between  1/10  and  1/7  of  the  natural  length 
scale . 

5 . 3  RESULTS . 

During  the  course  of  the  October  23  experiment,  the  burners 
were  run  for  ten  minutes,  with  a  "plume  steady  state  resulting 
after  five  minutes;  no  further  increase  of  the  top  plume  was 
visible"  [6].  For  comparison,  the  development  of  the  model  plume 
can  be  most  quickly  visualized  in  Fig.  19,  which  shows  the 
heights  reached  by  tracers  released  at  one  minute  intervals  from 
the  computational  cell  nearest  to  the  center  of  the  fire.  From 
this  figure  it  is  evident  that  the  model  plume  reached  its 
maximum  altitude  between  five  and  six  minutes  into  the  run,  in 
agreement  with  the  experiment.  Between  six  and  seven  minutes 


there  was  little  further  evolution  of  the  model  plume,  and  the 
fields  computed  at  seven  minutes  were  taken  to  represent  the 
plume  steady  state. 

The  top  plume  altitude  for  the  experiment  was  reported  to  be 
1260  meters,  while  the  maximum  altitude  of  the  plume  axis  was 
given  as  1050  meters  [6];  the  difference  between  these  two 
heights  is  due  to  the  inclination  of  the  axis  of  the  experimental 
plume  to  the  vertical,  and  may  be  regarded  as  a  rough  measure  of 
"experimental  error. "  The  two  reported  heights  have  been  plotted 
with  solid  horizontal  lines  in  Fig.  19,  where  it  can  be  seen  that 
the  model  plume  top  height  was  well  within  the  experimental 
range . 

In  both  the  experiment  and  the  model  simulation,  the  plume 
penetrated  the  inversion  and  rose  a  considerable  distance  into 
the  overlying  stable  layer.  This  process  is  depicted  in 
Figs.  39a-39f,  which  show  predicted  contours  of  the  model 
potential  temperature  field  at  one  minute  intervals,  starting 
from  the  initial  condition  depicted  in  Fig.  38a  (note  that  the 
highest  contour  plotted  in  these  figures  is  293°K) .  During  the 
first  two  minutes  the  plume  is  confined  to  the  lower  neutral 
layer,  although  disturbances  from  gravity-acoustic  waves  have 
propagated  into  the  overlying  stable  layer  (Figs.  39a-b) .  By 
three  minutes  into  the  run,  the  plume  has  reached  the  base  of  the 
inversion;  below  this  level,  entrainment  of  cooler  environmental 
air  causes  the  potential  temperature  of  the  plume  to  decrease 
with  height  (Fig.  39c).  As  the  plume  penetrates  into  the 
overlying  stable  layer  (Figs.  39d-e) ,  the  ambient  atmosphere  is 
displaced  by  a  large  region  of  uniform  potential  temperature, 
indicating  that  entrainment  of  environmental  air  towards  the  axis 
is  small  in  this  region.  The  air  in  the  plume  therefore  cools 
adiabatically  as  it  rises  into  the  stable  layer,  leading  to 
temperatures  near  the  axis  as  much  as  3.5°C  below  ambient  values 
at  an  altitude  of  1100  meters  (Fig.  39 f ) .  Temperatures  in  the 
overshoot  region  (above  750  m)  were  not  reported  in  [6],  so 


Cul'MNlt  *N0  TICHNCK.OC 


Model  potential  temperature  field  Figure  39b.  Model  potential  temperature  field 

at  1  minute  into  the  Meteotron  at  2  minutes  into  the  Meteotron 

simulation  run.  Contour  interval  simulation  run.  Contour  interval 


Figure  39e.  Model  potential  temperature  field  Figure  39f.  Model  potential  temperature  field 

at  5  minutes  into  the  Meteotron  at  6  minutes  into  the  Meteotron 

simulation  run.  Contour  interval  simulation  run.  Contour  interval 


this  feature  cannot  be  verified  directly.  During  airplane 
traversals  of  a  plume  generated  by  a  1000  MW  burn  at  the 
Meteotron  facility  on  July  8,  1978,  however,  temperatures  as  much 
as  2°C  below  ambient  were  detected  in  the  updraft  region  [12], 
lending  plausibility  to  the  under-temperatures  predicted  by  the 
model . 

Since  the  plume  was  cooler  than  the  ambient  atmosphere  above 
an  altitude  of  750  m,  its  penetration  into  the  upper  part  of  the 
stable  layer  is  due  to  inertia,  rather  than  to  buoyancy. 

Fig.  40a,  for  example,  shows  a  plot  of  the  velocity  field  at 
5  minutes,  when  the  plume  first  attains  its  maximum  altitude.  At 
this  time  positive  vertical  velocities  extend  up  to  an  altitude 
of  1200  m,  and  negative  vertical  velocities  are  virtually  absent. 
By  the  time  steady  state  is  achieved  at  seven  minutes  into  the 
run,  however,  (see  Figure  40b)  an  extensive  area  of  negative 
vertical  velocity  has  developed  as  the  overshooting  plume  falls 
back  towards  the  level  of  neutral  buoyancy.  Although  downdrafts 
near  the  upper  level  of  the  plume  were  not  reported  in  [6], 
airplane  traversals  of  several  experimental  plumes  conducted  at 
the  Meteotron  facility  during  the  summer  of  1978  did  detect 
downdrafts  of  several  meters/sec  [12],  lending  plausibility  to 
the  results  depicted  in  Fig.  40b. 

Fig.  20  shows  streamlines  of  the  steady-state  mass  flux 
field  generated  by  the  run.  The  flow  is  directed  parallel  to  the 
steamlines,  with  the  azimuthally  integrated  mass-flux  between  any 
two  neighboring  streamlines  equal  to  10 4  kg/sec;  note  that  below 
the  inversion,  horizontal  inflow  (entrainment)  causes  the 
vertical  mass  flux  of  the  plume  to  increase  with  height. 
Observations  of  the  visible  boundary  of  the  experimental  smoke 
plume  are  plotted  in  the  figure  with  black  circles.  While 
turbulent  diffusion  within  the  plume  causes  the  smoke  to  spread 
laterally,  horizontal  inflow  at  the  edges  of  the  plume  confines 
the  smoke  to  the  area  of  upward  motion  [6].  Fig.  20  shows  that 
good  agreement  between  the  edge  of  the  observed  smoke  plume  and 
the  boundary  of  the  simulated  inflow  area  has  been  obtained, 


indicating  that  the  width  of  the  plume  has  been  correctly 
predicted  by  the  model . 

Temperatures  within  the  plume  were  measured  with  a 
combination  of  kite  radiosonde  and  Barnes  radiometer.  The 
average  over-temperature  (i.e.,  excess  of  plume  temperature  over 
the  ambient  temperature  at  that  level)  for  eleven  Meteotron 
experiments  was  plotted  by  Benech  [6]  as  the  solid  line  in 
Fig.  21.  For  comparison,  model  over- temperatures  were  computed 
at  a  distance  from  the  axis  equal  to  one-half  the  observed  plume 
radius  at  the  level,  and  are  plotted  as  black  squares  in  Fig.  21. 
In  view  of  the  uncertainties  involved  in  selecting  a  single 
"representative"  plume  temperature  at  each  level,  satisfactory 
agreement  between  the  observed  and  computed  over-temperature  has 
been  obtained. 

5.4  DISCUSSION. 

From  the  above  comparisons  with  observed  data,  it  appears 
that  our  simulation  has  succeeded  in  capturing  the  essential 
dynamics  of  the  October  23  Meteotron  burn.  Since  the  aspect 
ratio  of  the  Meteotron  plume  was  considerably  larger  than  the 
aspect  ratios  of  the  firestorm  plumes  which  we  simulated  (compare 
Figs.  14  and  40),  physical  quantities  (e.g.,  mixing  lengths) 
cannot  be  scaled  directly  from  the  test-fire  to  firestorm  cases. 
In  spite  of  the  difference  in  aspect  ratios,  however,  a 
consistent  line  of  physical  reasoning  led  to  the  selection  of  an 
asymptotic  turbulent  length  scale  which  enabled  us  to  accurately 
calculate  the  stabilization  height  of  the  Meteotron  plume 
(see  Fig.  19).  Since  this  event  involved  the  penetration  of  a 
thermal  plume  into  an  overlying  stable  layer,  the  successful 
simulation  of  the  Meteotron  plume  encourages  us  to  believe  that 
the  DICE  code  is  capable  of  correctly  simulating  the  penetration 
of  a  firestorm  plume  into  the  stable  layers  of  the  stratosphere. 


SECTION  6 


THE  EFFECTS  OF  SWIRL 


6.1  OVERVIEW. 


Anecdotal  evidence  has  suggested  that  the  most  intense 
firestorms  created  during  World  War  II  were  accompanied  by 
strong,  large-scale  rotation,  and  some  investigators  have  claimed 
that  rotation  is  an  integral  part  of  the  firestorm  phenomenon, 
much  as  it  is  for  hurricanes  [7].  While  hurricane  flow  fields 
extend  for  hundreds  of  kilometers  and  draw  on  planetary  vorticity 
to  generate  their  rotation,  firestorm  flow  fields  extend  for  only 
10 's  of  kilometers,  and  must  rely  on  ambient  vorticity  in  the 
local  flow  field  (e.g.,  a  shear  in  the  horizontal  wind  speed)  to 
generate  large  rotational  velocities.  In  general,  rotational 
flow  (swirl)  may  affect  the  dynamics  of  thermal  plumes  either  by 
changing  the  mean  fields  (e.g.,  radial  pressure  gradients),  or  by 
changing  the  level  of  turbulence  (which  controls  the  entrainment 
of  ambient  air  into  the  plume) . 

6.2  EFFECTS  OF  SWIRL  ON  THE  MEAN  FLOW. 


In  order  to  investigate  the  effects  of  swirl  on  the  dynamics 
of  thermal  plumes,  the  DICE  code  has  been  extended  to  include  a 
tangential  velocity  component.  While  the  velocity  field  is  now 
"three-dimensional,"  axial  symmetry  has  been  preserved,  and  all 
quantities  are  still  functions  of  time,  radius,  and  height  only. 
Case  841  was  run  with  the  extended  version  of  the  code,  using  the 
the  same  parameters  and  turbulence  sub-model  as  Case  702  (see 
Table  2).  The  initial  surface  swirl  velocity  distribution  for 
this  case  was  given  by 

A 

w  =  2w  — ,  (19) 

l+xJ 
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where  x  =  x  /  10  km,  and  w  was  taken  to  be  10  m/sec.  The  profile 
of  the  initial  swirl  velocity  at  the  surface  is  shown  in  Fig.  22; 
in  the  free  atmosphere,  the  initial  swirl  velocity  decreases 
linearly  with  height  from  its  surface  value  to  zero  near  the 
tropopause . 

Figs.  23a-d  show  vector  plots  of  the  non-swirl  velocity 
component  and  contour  plots  of  the  swirl  speed  for  various  times 
in  the  run.  At  15  minutes  the  inflow  layer  is  not  yet  fully 
developed,  and  maximum  swirl  velocities  have  only  reached 
20  m/sec.  At  20  minutes  the  inflow  layer  has  almost  reached  the 
axis,  and  peak  swirl  speeds  near  the  axis  have  increased  to 
80  m/sec.  By  30  minutes  into  the  run  air  from  the  region  of 
maximum  initial  swirl  velocity  penetrates  to  the  axis,  and  the 
"spin-up"  of  these  velocities  has  generated  swirl  speeds  in 
excess  of  200  m/sec.  While  the  maximum  swirl  velocities  near  the 
surface  are  located  along  the  axis,  a  vortex  has  formed  above 
10  km,  with  maximum  swirl  speeds  located  about  1  km  from  the 
axis . 


A  time  history  of  the  maximum  swirl  velocity  in  the  grid 
(see  Fig.  41)  shows  that  the  "spin  up"  proceeds  quite  rapidly, 
once  air  from  the  edges  of  the  fire  (where  the  initial  swirl 
velocity  is  the  largest)  reaches  the  center  of  the  burning  region 
at  about  20  minutes  into  the  run.  At  later  times,  as  air  from 
regions  beyond  the  peak  initial  swirl  angular  momentum  reaches 
the  axis,  the  maximum  swirl  velocity  decreases,  reaching  160 
m/sec  by  45  minutes  into  the  run  (see  Fig.  23d);  by  this  time 
also  the  maximum  swirl  speeds  in  the  vortex  above  10  km  have 
collapsed  back  towards  the  axis. 

Since  this  run  did  not  include  the  effects  of  surface 
stress,  air  which  is  advected  towards  the  origin  in  the  inflow 
layer  will  tend  to  conserve  its  initial  angular  momentum  about 
the  axis  of  the  fire.  A  ring  of  air  moving  inwards  from 
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Figure  41.  Time  history  of  the  maximum  swirl  velocity  for  Case  841 


10  km  to  .33  km  (the  radial  distance  of  the  innermost  velocity 
grid-point),  for  example,  will  experience  an  increase  in  swirl 
velocity  from  a  value  of  10  m/sec  to  about  300  m/sec  if  no 
stresses  are  acting  on  it.  Stresses  exerted  on  the  air  in  the 
inflow  layer  by  air  above  it,  however,  limited  the  peak  swirl 
speeds  achieved  in  this  run  to  about  200  m/sec.  Of  course 
smaller  peak  swirl  speeds  would  be  obtained  if  surface  stresses 
had  been  included  in  the  model,  or  if  a  smaller  initial  swirl 
velocity  had  been  prescribed.  It  should  be  noted  that  even  for  a 
resting  initial  atmosphere,  the  Coriolis  effect  would  be  expected 
to  generate  swirl  velocities  about  an  order  of  magnitude  less 
than  those  obtained  in  this  run  (for  mid-latitude  sites). 

Theoretical  considerations  lead  us  to  believe  that  "in  a 
nonrotating  convective  cell,  the  vertical  velocity  would 
overshoot  its  equilibrium  altitude ...  in  contrast,  the  swirling 
updraft  will  have  an  adverse  pressure  gradient .. .which  not  only 
blocks  the  overshoot,  but  almost  certainly  forces  a  downdraft 
along  the  axis"  ([9],  p.  138).  For  our  baseline  non-rotating 
case  (702),  the  plume  shows  considerable  overshoot  by  20  minutes 
(see  Fig.  14b);  for  the  rotating  case  (841),  a  plot  of  the  non¬ 
swirl  velocity  component  at  20  minutes  (see  Fig.  23b)  showB  that 
the  rise  of  the  plume  has  been  substantially  blocked,  and  that  an 
area  of  negative  vertical  velocity  has  indeed  formed  along  the 
axis.  Although  the  swirling  plume  does  penetrate  into  the 
stratosphere  at  later  times  (see  Figs.  23c-d),  the  plume  heights 
remain  about  3  km  lower  than  the  corresponding  heights  for  Case 
702  (compare  with  Figs.  14c-d).  Due  to  the  reduced  overshoot, 
the  tracer  cloud  for  case  841  (see  Fig.  42)  does  not  show  the 
secondary  region  of  lofted  soot  which  was  obtained  in  case  702  at 
higher  altitudes  near  the  axis  (see  Fig.  15). 

The  results  of  Case  841  show  that  large  rotational 
velocities  may  develop  in  the  interior  of  a  10-km  radius  fire, 
given  an  initial  ambient  swirl  of  10  m/sec  at  the  fire  edge,  and 
that  the  effect  of  the  swirl  is  to  delay  plume  rise  and  lower  the 
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stabilization  height  somewhat.  Interaction  between  the  swirl 
velocity  field  and  the  turbulence  was  not  considered  in  this 
case,  however,  and  it  is  possible  that  the  swirl  may  act  to 
suppress  turbulence,  much  as  a  stable  density  stratification 
does.  As  was  shown  by  the  comparison  b  tween  Cases  702  and  705 
in  Section  4,  this  would  tend  to  intensify  the  plume. 

6.3  EFFECTS  OF  SWIRL  ON  TURBULENCE. 

A  given  distribution  of  swirl  velocity  may  act  as  a  source 
or  sink  of  turbulent  kinetic  energy  ( TKE ) .  When  only  swirl 
velocity  is  present,  a  local  balance  between  production  and 
dissipation  terms,  similar  to  that  used  in  the  turbulence  sub¬ 
model  described  in  Section  4,  may  be  solved  to  give  a  diagnostic 
relation  for  the  TKE  in  terms  of  the  mean  swirl  velocity 
distribution  [10].  This  relation  turns  out  to  have  a  form 
similar  to  that  which  would  be  given  by  a  simple  mixing  length 
argument;  the  resulting  turbulent  stresses,  however,  also  depend 
on  the  curvature  Richardson  number 


c  0_ 
8r 


through  the  stability  function 


(l-Rc) 


Here  w  is  the  swirl  velocity,  r  the  radial  distance,  and  K  a 
fixed  constant. 


A  plot  of  the  stability  function  SC(RC),  taken  from  [10],  is 
shown  in  Fig.  24.  For  values  of  Sc  greater  than  1,  the  curvature 
effect  amplifies  the  turbulence,  while  for  values  of  Sc  less  than 
1,  the  curvature  has  a  damping  effect.  For  values  of  Rc  between 


0.185  and  5.416,  Sc  becomes  negative;  in  this  case  destruction  of 
TKE  by  the  mean  swirl  velocity  field  exceeds  production,  and  no 
turbulence  exists. 

When  the  generation  of  TKE  by  non-swirl  winds  and  buoyancy 
is  taken  into  account,  the  destruction  of  TKE  by  the  swirl 
velocity  field  may  still  exert  a  damping  effect  on  the  turbulence 
or  extinguish  it  completely,  depending  on  the  relative  magnitudes 
of  the  source  and  sink  terms  in  the  TKE  budget  equation. 
Unfortunately  this  equation  is  too  complicated  in  the  general 
case  to  be  solved  diagnostically  for  the  TKE,  as  was  done  for  the 
case  with  only  non-swirl  terms  (Section  4),  or  for  the  case  with 
only  swirl  terms  [10].  To  adequately  treat  the  case  where  swirl 
and  non-swirl  terms  are  active  simultaneously,  it  becomes 
necessary  to  include  separate  predictive  equations  in  the  model 
for  all  components  of  the  turbulent  stress  tensor,  as  well  as  for 
the  turbulent  heat  fluxes. 

While  no  such  scheme  has  yet  been  implemented  in  the  DICE 
code,  a  model  developed  by  the  Aeronautical  Research  Associates 
of  Princeton,  Inc.  (ARAP)  currently  has  this  capability.  In  a 
recent  study  of  tornado  dynamics  [9],  two  versions  of  the  ARAP 
code  were  used;  one  which  had  a  simplified  turbulence  model 
similar  to  that  used  in  DICE  Cases  702  and  705  (which  does  not 
include  the  effects  of  curvature  on  turbulence),  and  one  which 
used  the  fully  predictive  model  described  above  (which  implicitly 
includes  the  effects  of  curvature  on  turbulence) .  Non- 
dimensional  plots  of  the  turbulent  eddy  viscosity  field  produced 
by  the  runs  with  the  simplified  and  full  turbulence  models  are 
shown  in  Figs.  25a  and  25b,  respectively.  Near  the  surface  the 
turbulence  profiles  are  similar,  but  at  higher  altitudes  the  full 
turbulence  model  produces  much  less  eddy  viscosity,  indicating 
that  rotational  flow  can  have  a  damping  effect  on  turbulence 
above  the  boundary  layer. 
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Since  the  ARAP  study  dealt  with  swirling  plumes,  the  results 
are  of  interest  to  the  rotating  firestorm  case,  although  the 
tornado  problem  includes  no  buoyant  generation  of  TKE  in  the 
boundary  layer.  Some  preliminary  calculations  done  at  CRT,  in 
fact,  indicate  that  buoyant  production  of  TKE  in  the  boundary 
layer  of  a  firestorm  would  proceed  at  a  considerably  higher  rate 
than  its  destruction  by  a  plausible  swirl  velocity  distribution. 
Together  with  the  result  obtained  in  the  ARAP  study,  this 
suggests  that  the  suppression  of  turbulence  by  swirling  flow  in 
the  boundary  layer  of  a  firestorm  would  probably  be  small. 

In  comparing  Cases  702  (£»  =  100  m)  and  705  (£»  =  25  m)  in 
Section  4,  we  saw  that  reducing  the  turbulence  by  decreasing  the 
mixing  length  from  100  m  to  25  m  throughout  the  computational 
domain  resulted  in  a  large  increase  in  the  plume  stabilization 
height.  While  it  appears  unlikely  that  swirl  would  decrease  the 
turbulence  in  the  boundary  layer,  the  ARAP  results  suggest  that 
swirl  could  significantly  reduce  the  level  of  turbulence  in  the 
region  above  the  boundary  layer.  In  order  to  estimate  the  likely 
effects  of  swirl  on  the  plume,  then,  we  ran  Case  846  with  the 
same  parameters  as  Case  702  (see  Table  2),  keeping  the  mixing 
length  fixed  at  100  m  between  the  surface  and  0.5  km  altitude, 
but  letting  it  decrease  linearly  to  25  m  at  1  km  altitude  (see 
Fig .  26 ) . 

Figs.  27a-b  show  the  resulting  velocity  fields  at  20  and 
30  minutes.  Comparing  with  the  corresponding  plots  for  Case  702 
(see  Figs.  14b-c),  we  see  that  cutting  back  the  mixing  length 
from  100  m  to  25  m  above  the  boundary  layer  causes  an  enhancement 
in  the  plume's  height  of  about  3  km.  A  much  larger  enhancement, 
however,  is  gained  by  reducing  the  mixing  length  to  25  m 
throughout  the  computational  domain,  as  was  shown  by  the 
comparison  between  Cases  702  and  705.  This  result  is  not 
unexpected,  in  view  of  the  conclusion  reached  in  Section  4  that 
entrainment  within  the  boundary  layer  determines  the  core 


temperature,  and  ultimately  the  stabilization  height,  of  the 
plumes  from  large-area  fires. 

6.4  DISCUSSION. 

Case  841,  described  in  Section  6-2,  showed  that  high  swirl 
velocities  can  develop  near  the  core  of  a  large-area  fire,  given 
a  moderate  initial  ambient  swirl  velocity  field.  We  saw  that  the 
large  swirl  velocities  which  developed  were  due  to  the 
approximate  conservation  of  angular  momentum  in  the  inflow  layer, 
which  can  be  expected  to  hold  as  long  as  the  flow  field  maintains 
axial  symmetry.  By  comparing  with  a  similar  non-swirl  case,  we 
also  saw  that  the  effects  of  swirl  on  the  mean  flow  dynamics  were 
to  delay  plume  rise  somewhat,  and  to  lower  the  plume  height  by 
about  3  km. 

The  question  of  how  swirl  will  affect  the  turbulence  in  a 
large-area  fire  is  much  more  difficult,  and  at  present  we  can 
only  make  reasonable  estimates  of  this  effect.  By  assuming  that 
swirl  decreases  the  turbulence  in  the  region  above  the  boundary 
layer,  we  obtained  a  moderate  enhancement  in  the  plume  height. 
Suppression  of  turbulence  in  the  inflow  layer  cannot  be  ruled 
out,  however,  particularly  in  view  of  the  large  rotational 
velocities  which  develop  near  the  axis.  The  definitive 
resolution  of  this  question  awaits  the  extension  of  the  DICE  code 
to  predict  the  full  set  of  turbulent  fluxes,  or  the  application 
of  a  model  which  already  has  that  capability,  such  as  the  ARAP 
code,  to  the  firestorm  problem. 
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